text:
# PHS3350
# Week 5 - Energy levels of a family of non-Hermitian Hamiltonians
# "what I cannot create I do not understand" - R. Feynman.
# Ana Fabela Hinojosa, 04/04/2021

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import fsolve
from scipy.special import gamma
import scipy.special as sc
from scipy import linalg
from tqdm import tqdm
from odhobs import psi as cpsi_blank

plt.rcParams['figure.dpi'] = 200
np.set_printoptions(linewidth=200)

E0 = 1.1563

def complex_quad(func, a, b, **kwargs):
    # Integration using scipy.integratequad() for a complex function
    def real_func(*args):
        return np.real(func(*args))

    def imag_func(*args):
        return np.imag(func(*args))

    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return real_integral[0] + 1j * imag_integral[0]

## TEST
def complex_fsolve(func, E0, **kwargs):
    # root finding algorithm. FINDS: Energy values from error() function
    # call: complex_fsolve(error, E0, args=(ϵ, n))
    def real_func(*args):
        return np.real(func(*args))

    real_root = fsolve(real_func, E0, **kwargs)
    #     def imag_func(*args):
    #         return np.imag(func(*args))
    #     imag_root = fsolve(imag_func, E0, **kwargs)
    # Check that the imaginary part of func is also zero
    value = func(real_root[0], *kwargs['args'])
    assert abs(np.imag(value)) < 1e-10, "Imaginary part wasn't zero"
    # print(f"E = {Energies[0]:.04f}")
    return real_root[0]

def integrand(x_prime, tp_minus, E, ϵ):
    # Change of variables integrand
    x = x_prime + 1j * np.imag(tp_minus)
    return np.sqrt(E - x ** 2 * (1j * x) ** ϵ)

def LHS(n):
    # Quantization condition
    return (n + 1 / 2) * np.pi

def RHS(E, ϵ):
    # Integral defining E
    # integration LIMITS
    tp_minus = E ** (1 / (ϵ + 2)) * np.exp(1j * np.pi * (3 / 2 - (1 / (ϵ + 2))))
    tp_plus = E ** (1 / (ϵ + 2)) * np.exp(-1j * np.pi * (1 / 2 - (1 / (ϵ + 2))))
    tp_minus_prime = np.real(
        E ** (1 / (ϵ + 2)) * np.exp(1j * np.pi * (3 / 2 - (1 / (ϵ + 2))))
        - 1j * np.imag(tp_minus)
    )
    tp_plus_prime = np.real(
        E ** (1 / (ϵ + 2)) * np.exp(-1j * np.pi * (1 / 2 - (1 / (ϵ + 2))))
        - 1j * np.imag(tp_minus)
    )
    # print(tp_minus_prime)
    # print(tp_plus_prime)
    return complex_quad(integrand, tp_minus_prime, tp_plus_prime, args=(tp_minus, E, ϵ))

def error(E, ϵ, n):
    return RHS(E, ϵ) - LHS(n)

def compare():
    # comparison to WKB results reported for  (ϵ, n) = (1, 0) using IC: E0 = E0
    E_RK = E0
    E_WKB = 1.0943
    diff_RK_WKB = E_RK - E_WKB
    diff_RK_mine = E_RK - complex_fsolve(error, E0, args=(1, 0))
    how_many_sigmas_theory = diff_RK_WKB / E_RK
    how_many_sigmas_mine = diff_RK_mine / E_RK
    print(
        f"\nnumber of 𝞼 away is WKB from exact result: {abs(how_many_sigmas_theory):.3f}"
    )
    print(f"number of 𝞼 away am I from exact result: {abs(how_many_sigmas_mine):.3f}\n")

def analytic_E(ϵ, n):
    # Bender equation (34) pg. 960
    top = gamma(3 / 2 + 1 / (ϵ + 2)) * np.sqrt(np.pi) * (n + 1 / 2)
    bottom = np.sin(np.pi / (ϵ + 2)) * gamma(1 + 1 / (ϵ + 2))
    return (top / bottom) ** ((2 * ϵ + 4) / (ϵ + 4))

def brute_force(func, E, ϵ):
    # BRUTE INTEGRAL
    def real_func(*args):
        return np.real(func(*args))

    tp_minus = E ** (1 / (ϵ + 2)) * np.exp(1j * np.pi * (3 / 2 - (1 / (ϵ + 2))))
    tp_plus = E ** (1 / (ϵ + 2)) * np.exp(-1j * np.pi * (1 / 2 - (1 / (ϵ + 2))))
    tp_minus_prime = E ** (1 / (ϵ + 2)) * np.exp(
        1j * np.pi * (3 / 2 - (1 / (ϵ + 2)))
    ) - 1j * np.imag(tp_minus)
    tp_plus_prime = E ** (1 / (ϵ + 2)) * np.exp(
        -1j * np.pi * (1 / 2 - (1 / (ϵ + 2)))
    ) - 1j * np.imag(tp_minus)
    # domain & differential (infinitesimal)
    x_prime = np.linspace(tp_minus_prime, tp_plus_prime, 50000)
    dx_prime = x_prime[1] - x_prime[0]
    return np.sum(real_func(x_prime, E, ϵ) * dx_prime)

## Runge-Kutta, finding IC!
def find_k(x, ϵ, E):
    return np.sqrt(x ** 2 * (1j * x) ** ϵ - E)

# Schrödinger equation
def Schrodinger_eqn(x, Ψ):
    psi, psi_prime = Ψ
    psi_primeprime = (x ** 2 * (1j * x) ** ϵ - E) * psi
    Ψ_prime = np.array([psi_prime, psi_primeprime])
    return Ψ_prime

def Runge_Kutta(x, delta_x, Ψ):
    k1 = Schrodinger_eqn(x, Ψ)
    k2 = Schrodinger_eqn(x + delta_x / 2, Ψ + k1 * delta_x / 2)
    k3 = Schrodinger_eqn(x + delta_x / 2, Ψ + k2 * delta_x / 2)
    k4 = Schrodinger_eqn(x + delta_x, Ψ + k3 * delta_x)
    return Ψ + (delta_x / 6) * (k1 + 2 * k2 + 2 * k3 + k4)

#######################################function calls############################################

# WKB approximation
# IC based on RK results given (ϵ, n) = (1, 0)
ϵ = 1
E1 = 1.1563
E = E1

N = 100
epsilons = np.linspace(-1.0, 0, N)

tp_minus = E ** (1 / (ϵ + 2)) * np.exp(1j * np.pi * (3 / 2 - (1 / (ϵ + 2))))
tp_plus = E ** (1 / (ϵ + 2)) * np.exp(-1j * np.pi * (1 / 2 - (1 / (ϵ + 2))))

tp_minus_prime = E ** (1 / (ϵ + 2)) * np.exp(
    1j * np.pi * (3 / 2 - (1 / (ϵ + 2)))
) - 1j * np.imag(tp_minus)
tp_plus_prime = E ** (1 / (ϵ + 2)) * np.exp(
    -1j * np.pi * (1 / 2 - (1 / (ϵ + 2)))
) - 1j * np.imag(tp_minus)

######################## FINAL plot of eigenvalues ###################################

# Energies_2 = []
# for n in range(10):
#     E_s = []
#     for ϵ in np.linspace(0, 3, 30):
#         E = complex_fsolve(error, E1, args=(ϵ, n))
#         E_s.append(E)
#         ## print(f" {ϵ = }, {n = }, {E = }")
#     Energies_2.append(E_s)
# np.save("Energies_unbroken.npy", Energies_2)
Energies_2 = np.load("Energies_unbroken.npy")

def figure_final_form():
    for E_ϵs in Energies_2:
        # print(E_ϵs)
        ϵ = np.linspace(0, 3, 30)
        plt.plot(ϵ, E_ϵs, "o-", color='k', markersize=1) # MARK OR UNMARK

    eigenvectors_list = []
    for i, ϵ in enumerate(epsilons):
        matrix = np.load(f'matrices/matrix_{i:03d}.npy')
        eigenvalues, eigenvectors = linalg.eig(matrix)
        eigenvectors_list.append(eigenvectors)

        # full figure # MARK OR UNMARK
        positive_evals = [
            i for i in eigenvalues if 0 < np.real(i) < 20 and abs(np.imag(i)) < 0.3
        ]

        # # # broken region # MARK OR UNMARK
        # positive_evals = [
        #     i for i in eigenvalues if 0 < np.real(i) < 20 and abs(np.imag(i)) < 6
        # ]

        sorted_eigenvalues = sorted(positive_evals, key=lambda x: np.real(x))
        sorted_eigenvalues = sorted_eigenvalues[:11]

        # ϵ_list = np.full(len(sorted_eigenvalues), ϵ)

        # mask_imag = 1e-6 < abs(np.imag(sorted_eigenvalues)) # MARK OR UNMARK
        # plt.plot(
        #     ϵ_list[mask_imag],
        #     np.imag(sorted_eigenvalues)[mask_imag],
        #     marker='.',
        #     linestyle='None',
        #     color='r',
        #     markersize=1,
        # )

        # plt.plot(
        #     ϵ_list,
        #     np.real(sorted_eigenvalues),
        #     marker='.',
        #     linestyle='None',
        #     color='k',
        #     markersize=2.5,
        # )

        # mask_real = 1e-6 < abs(np.imag(sorted_eigenvalues))# MARK OR UNMARK
        # plt.plot(
        #     ϵ_list[mask_real],
        #     np.real(sorted_eigenvalues)[mask_real],
        #     marker='.',
        #     linestyle='None',
        #     color='xkcd:azure',
        #     markersize=2,
        # )

    # # full figure # MARK OR UNMARK
    # plt.axis(xmin=-1, xmax=3, ymin=0, ymax=20)
    # plt.axvline(0, color='limegreen', linestyle=':', label="PT-symmetry breaking")
    # plt.legend()

    # # # only broken symmetry region # MARK OR UNMARK
    # # plt.axis(xmin=-1, xmax=0, ymin=-2, ymax=12)
    # # plt.axhline(0, color='grey', linestyle='-')

    # plt.xlabel("ϵ")
    # plt.ylabel("E")
    # plt.savefig("NHH_eigenvalues.png") # MARK OR UNMARK
    # plt.show()
    return np.array(eigenvectors_list)

coefficients = figure_final_form()

# ####################### Eigenvectors plot ##################################

# def spatial_wavefunctions(N, x, epsilons):
#     x[x == 0] = 1e-200
#     PSI_ns = []
#     for n in range(N):
#         psi_n = cpsi_blank(n, x)
#         PSI_ns.append(psi_n)
#     PSI_ns = np.array(PSI_ns)
#     np.save(f"PSI_ns.npy", PSI_ns)

#     eigenstates = []
#     for i, ϵ in enumerate(epsilons):
#         c = coefficients[i]
#         for j in range(N):  # for each eigenvector
#             d = c[:, j]
#             psi_jx = np.zeros(x.shape, complex)
#             for n in range(N):  # for each H.O. basis vector
#                 psi_jx += d[n] * PSI_ns[n]
#             plt.plot(x, abs(psi_jx) ** 2)
#             plt.savefig(f"spatial_wavefunctions/wavefunction_{ϵ}_{n:03d}.png")
#             plt.clf()
#             eigenstates.append(psi_jx)

#     np.save(f'eigenstates.npy', np.array(eigenstates))

# N = 100
# epsilons = np.linspace(-1.0, 0, N)
# xs = np.linspace(-20, 20, 1024)
# spatial_wavefunctions(N, xs, epsilons)

# ################# another one ###############################################
def spatial_wavefunctions_orthogonality_check(N, x, epsilons):
    x[x == 0] = 1e-200
    delta_x = x[1]-x[0]
    PSI_ns = np.load("PSI_ns.npy")

    eigenstates = []
    for i, ϵ in enumerate(epsilons):
        if i != 30:
            continue
        # if i != 60:
        #     continue
        c = coefficients[i]
        for j in range(N):  # for each eigenvector
            d = c[:, j]
            psi_jx = np.zeros(x.shape, complex)

            for n in range(N):  # for each H.O. basis vector
                psi_jx += d[n] * PSI_ns[n]

            if j == 39:
                psi_39 = psi_jx
            if j == 40:
                psi_40 = psi_jx

            plt.plot(x, abs(psi_jx) ** 2)
            print(f"saving {ϵ = }, {j = }")
            plt.savefig(f"some_wavefunctions_-0.7/wavefunction_{ϵ}_{j:03d}.png")
            # plt.savefig(f"some_wavefunctions_-0.4/wavefunction_{ϵ}_{j:03d}.png")
            plt.clf()
            eigenstates.append(psi_jx)
        print(np.vdot(psi_39, psi_40) * delta_x)
    np.save(f'eigenstates-0.7.npy', np.array(eigenstates))
    # np.save(f'eigenstates-0.4.npy', np.array(eigenstates))

plt.clf()
N = 100
xs = np.linspace(-20, 20, 2048)
spatial_wavefunctions_orthogonality_check(N, xs, epsilons)

####################### Eigenvectors plot ##################################

###################### Runge-Kutta test call ###############################

# psi = np.empty_like(xs, dtype = "complex_")
# for i in range(len(xs)):
#     psi[i] = Runge_Kutta(x, delta_x, Ψ0)[0]
#     # print(psi[i])
#     x += delta_x

##################### Runge-Kutta test call ###############################

# ######################### WKB TEST 1 ####################################
# def whats_up_with_integrand(x_values, E, ϵ):
#     # checking the integration path of the integrand in the x-complex plane
#     reals = []
#     imaginary = []
#     for x in x_values:
#         complex_num = np.sqrt(E - x**2 * (1j * x)**ϵ)

#         reals.append(np.real(complex_num))
#         imaginary.append(np.imag(complex_num))

#     plt.plot(reals, imaginary, '-')
#     plt.plot(reals[0], imaginary[0],'go', label='start here')
#     plt.plot(reals[5], imaginary[5],'-o', markersize='1.2')

#     plt.plot(reals[-1], imaginary[-1],'ro', markersize='1.2', label='finish here')
#     plt.legend()
#     plt.ylabel(r'$Im(\sqrt{E - x^2 (i x)^\epsilon})$')
#     plt.xlabel(r'$Re(\sqrt{E - x^2 (i x)^\epsilon})$')

#     # plt.title("Bender's integration contour")
#     plt.title("change or variables integration contour")
#     plt.show()

# # Bender's integral
# x_values = np.linspace(tp_minus, tp_plus, 10000)

# # change of variables integral
# x_values = np.linspace(tp_minus_prime, tp_plus_prime, 10000) + 1j * np.imag(tp_minus)

# whats_up_with_integrand(x_values, E0, ϵ)
# ######################### WKB TEST 1 #####################################

########################### WKB TEST 2 #####################################
# # ITERATIVE approach 1 for ϵ = 1
# Energies_1 = []
# for n in range(10):
#   E = complex_fsolve(error, E1, args=(1, n))
#   Energies_1.append(E)

# # comparison to WKB and RK reported in Bender
# n = range(10)
# E_RK= [1.1563, 4.1093, 7.5623, 11.3144, 15.2916, 19.4515, 23.7667, 28.2175, 32.7891, 37.4698]
# E_WKB = [1.0943, 4.0895, 7.5489, 11.3043, 15.2832, 19.4444, 23.7603, 28.2120, 32.7841, 37.4653]
# plt.plot(n, Energies_1 , label="my calculated energies")
# plt.plot(n, E_RK , label=r"$E_{RK}$")
# plt.plot(n, E_WKB , label=r"$E_{WKB}$")
# plt.legend()
# plt.xlabel("n")
# plt.ylabel("E")
# plt.show()
########################## WKB TEST 2 #####################################

######################## WKB unbroken region ##############################

# # ITERATIVE
# Energies = []
# for n in range(10):
#     E_ϵs = []
#     for ϵ in np.linspace(0, 3, 30):
#         E_ϵ = complex_fsolve(error, E0, args=(ϵ, n))
#         E_ϵs.append(E_ϵ)
#     Energies.append(E_ϵs)

# #PLOTING ITERATIVE
# for E_ϵs in Energies:
#     ϵ = np.linspace(0, 3, 30)
#     plt.plot(ϵ, E_ϵs, "o-", markersize=2)
# plt.ylim(0, 20)
# plt.xlabel("ϵ")
# plt.ylabel("E")
# plt.show()

######################## WKB unbroken region ##############################

max_stars_repo_name: anabanami/Bender

max_issues_repo_name: anabanami/Bender

max_forks_repo_name: anabanami/Bender

reason: import numpy,import scipy,from scipy


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

# # PRÁCTICA 2: tipologia y ciclo de vida de los datos

# LAURA BELENGUER QUEROL Y LEIRE ALEGRIA MURILLO

# In[1]:

#Importing required packages
import numpy as np
import pandas as pd

### Graphic libraries
import matplotlib.pyplot as plt
import seaborn as sns

import warnings
warnings.filterwarnings('ignore')

# ## 1. Descripción del dataset

# Leemos los datos de la librería Kaggle indicada en el ejercicio
# Tenemos 1599 observaciones y 12 variables

# In[2]:

# Cargamos los datos
wine = pd.read_csv('winequality-red.csv')

# In[3]:

# Dimension del dataset
wine.shape

# In[4]:

# Observamos nuestros datos y comprobamos que esten bien cargados
wine.head()

# ### 1.1. Análisis de la variable target

# Su valor está entre 0 y 10. Observamos que solo tenemos valores entre 3 y 8. La mayoría de los vinos tienen una calificacion de 5 o 6.

# In[5]:

wine['quality'].unique()

# In[6]:

# Target distribution
wine['quality'].hist(bins=40, edgecolor='b', linewidth=1.0,
              xlabelsize=8, ylabelsize=8, grid=False, figsize=(6,2), color='red')
plt.tight_layout(rect=(0, 0, 1.2, 1.2))
plt.suptitle('Red Wine Quality Plot', x=0.65, y=1.25, fontsize=14);

# Si un vino es excelente tiene una nota de 7 o superior, por lo que mas adelante creamos una variable llamada 'rating' donde clasificamos los vinos entre 'superior','inferior' o 'correcto'. Si su nota es superior a 7 sera superior, de lo contrario sera inferior.

# In[7]:

wine.quality.value_counts().sort_index()

# In[8]:

sns.countplot(x='quality', data=wine)

# ### 1.2 Análisis del resto de variables

# Mostramos los valores estadísticos de cada atributo.

# In[9]:

#Informacion sobre las columnas
wine.describe()

# Calculámos el histograma de cada atributo.

# In[10]:

target_wine = wine['quality']
data_wine = wine.drop(columns = ['quality'])

# In[11]:

# Histogramas de todas las variables
data_wine.hist(bins=40, edgecolor='b', linewidth=1.0,
                          xlabelsize=8, ylabelsize=8, grid=False,
                          figsize=(16,6), color='red')
plt.tight_layout(rect=(0, 0, 1.2, 1.2))
plt.suptitle('Red Wine Univariate Plots', x=0.65, y=1.25, fontsize=14);

# ## 2. Integración y selección de los datos de interés a analizar

# ### 2.1. Integración

# En este caso los datos se encuentran todos en el mismo dataset, por lo tanto no hay que integrar datos de otro dataset.

# Renombramos los nombres de los atributos para eliminar los espacio en blanco

# In[12]:

wine.rename(columns={'fixed acidity': 'fixed_acidity','citric acid':'citric_acid','volatile acidity':'volatile_acidity','residual sugar':'residual_sugar','free sulfur dioxide':'free_sulfur_dioxide','total sulfur dioxide':'total_sulfur_dioxide'}, inplace=True)
wine.head(n=5)

# ### 2.2. Selección

# Una de las primeras etapas en el preprocesado de los datos es el filtrado o selección de datos de interés.
#
# En esta fase también es habitual realizar una exploración de los datos.
#
# El preprocesado de los datos también puede incluir la creación de nuevas variables a partir de la extracción de características de los datos originales. Realizaremos está reducción con la variable 'quality'.

# In[13]:

# Basado en las categorías encontradas en el apartado anterior para 'quality'
# creamos una nueva variable 'rating' con 3 categorias para los vinos
conditions = [
    (wine['quality'] > 5),
    (wine['quality'] <= 5)
]
rating = ['1', '0']
wine['rating'] = np.select(conditions,rating)
wine.rating.value_counts()

# In[14]:

wine.groupby('rating').mean()

# ## 3. Limpieza de los datos

# ### 3.1. ¿Los datos contienen ceros o elementos vacíos? ¿Cómo gestionarías cada uno de estos casos?

# In[15]:

# Comprobamos si faltan valores "missing values"
print(wine.isna().sum())

# In[16]:

# Comprobamos si hoy valores nulos
print(wine.isnull().sum())

# En caso que hubiesen datos NaN en alguna de las columnas, utilizariamos: wine.dropna() para eliminar aquellas filas que contienen dichos valores NaN.

# ### 3.2. Identificación y tratamiento de valores extremos

# Los valores extremos (extreme scores o outliers) son aquellos datos que se encuentran muy alejados de la distribución normal de una variable o población.
#
# Son observaciones que se desvían tanto del resto que levantan sospechas sobre si fueron generadas mediante el mismo mecanismo.
#
# Estos valores pueden afectar de forma adversa los resultados de los análisis posteriores, al incrementar el error en la varianza de los datos y sesgar significativamente los cálculos y estimaciones.

# En el apartado 1.2 hemos visto que algunas variables/atributos tienen una desviación standard (std) bastante elevada si observamos el valor de la media.
#
# - citric acid
#
# - volatile acidity
#
# - residual sugar
#
# - chlorides
#
# - free sulfur dioxide
#
# - total sulfur dioxide
#
# Con los gráficos boxplot analizamos las variables de manera visual.
#
# También basados en los valores de calidad, visualizamos como dichos parámetros afectan a la calidad del vino.

# In[17]:

# Mostramos el boxplot de las variables (menos 'rating')
plt.figure(figsize=(12,6))
ax = sns.boxplot(data=wine.drop(columns = ['quality','rating']), orient="h")

# In[18]:

#Analizamos las variables una por una
fig, axes = plt.subplots(nrows=2, ncols=6, figsize=(12,5))

for i, ax in enumerate(axes.flatten()):
    sns.boxplot(y= wine.iloc[:, i],  orient='v' , ax=ax)

fig.tight_layout()

# Observamos que 'casi todas las variables presentan valores extremos. Podríamos identificar como outliers aquellos puntos que se encuentran, unidimensionalmente a más de 2 desviaciones estándar de la media (Distancia Mahalanobis)

# In[19]:

data_wine_nout = data_wine[data_wine.apply(lambda x: np.abs(x - x.mean()) / x.std() < 2).all(axis=1)]
print(data_wine_nout.shape)

# Estadísticas sin outliers y sin variable quality/rating
data_wine_nout.describe()

# Verificamos los resultados comparándolos con los originales

# In[20]:

# Estadísticas datos originales sin la variable quality/rating
data_wine.describe()

# In[21]:

# Mostramos el boxplot de las variables (menos 'quality'/'rating')
plt.figure(figsize=(12,6))
ax = sns.boxplot(data=data_wine_nout, orient="h")

# Aplicando el criterio de Mahalanobis para eliminar los valores extremos hemos reducido el tamaño de nuestro dataset de 1599 muestras a 1124. Con los boxplots vemos que se han reducido considerablemente los outliers para la mayoría de variables. Aunque ciertas variables como 'residual sugar' y 'chlorides' presentan un alto número de valores por encima del tercer cuartil (75% de los datos).

# De todos modos eliminar dichos valores extremos puede dar lugar a eliminar aquellos vinos que marcan la diferencia por atributos excepcionales, por ello nos quedaremos con los atributos originales sin eliminar ningún outlier.

# #### Valores duplicados

# In[22]:

extra = wine[wine.duplicated()]
extra.shape

# Hay 240 valores duplicados, pero no los eliminamos ya que pueden deberse a calificaciones idénticas dadas por diferentes
# catadores.

# # 4. Análisis de los datos.

# ## 4.1. Selección de los grupos de datos que se quieren analizar/comparar (planificación de los análisis a aplicar).

# En los siguientes apartados realizaremos el análisis inferencial que tiene por objetivo modelar los datos a través de una distribución conocida. Partiendo de la premisa que el conjunto de datos estudiado representa una fracción de la totalidad de una población, su objetivo es inferir cómo es esa población, asumiendo un grado de error en las estimaciones por el hecho de disponer de una muestra reducida de los datos.

# ## 4.2. Comprobación de la normalidad y homogeneidad de la varianza.

# ***Normalidad***

# Para la comprobación de que los valores que toman nuestras variables cuantitativas provienen de una población distribuida normalmente, utilizaremos la prueba de normalidad de Shapiro.
#
# Se comprueba que para que cada prueba se obtiene un p-valor superior al nivel de signiﬁcación preﬁjado α = 0.05. Si esto se cumple, entonces se considera que dicha variable sigue una distribución normal.

# H_0 = Tiene una distribución normal
#
# H_1 = No tiene una distribución normal

# p <= alpha: rechazamos H_0, no hay normalidad.
#
# p > alpha: no rechazamos H_0, hay normalidad.

# In[23]:

from numpy.random import seed
from numpy.random import randn

# In[24]:

from scipy.stats import shapiro

# Los atributos y target
names = list(wine.columns)

# no consideramos 'rating'
for i in range(len(names)-1):
    stat, p = shapiro(wine.iloc[:, i])
    print(names[i], ': p = %.3f' % (p))

# Todas las variables tienen el p valor inferior a 0.05, por lo tanto, ninguna variable sigue una distribución normal.

# ***Homogeneidad***

# Cuando los datos no cumplen con la condición de normalidad, cómo hemos visto en el apartado anterior, utilizamos el test Fligner-Killeen
#
# La hipótesis nula asume igualdad de varianzas en los diferentes grupos de datos, por lo que p-valores inferiores al nivel de significancia (0.05) indicán heterocedasticidad.
#
# En este caso lo aplicaremos a todas las variables, comparando los atributos (data_wine) con la variable target (quality)

# In[25]:

# Aplicamos la prueba de Fligner
from scipy.stats import fligner

names = list(data_wine.columns)

# Consideramos las variables/atributos y la comparamos con quality
for i in range(len(names)):
    stat, p = fligner(data_wine.iloc[:, i], wine['quality'])
    print(names[i], ': p = %.3f' % (p))

# Encontramos por tanto Heterocedastidad

# ## 4.3. Aplicación de pruebas estadísticas para comparar los grupos de datos.
# #### En función de los datos y el objetivo del estudio, aplicar pruebas de contraste de hipótesis, correlaciones, regresiones, etc. Aplicar al menos tres métodos de análisis diferentes.

# ### ***Contraste de hipótesis***

# * #### **Test Mann Whitney**

# En el caso que consideremos 2 niveles para la variable 'rating' o podemos utilizar también la variable 'quality' directamente

# In[26]:

from scipy.stats import mannwhitneyu

for i in range(len(names)):
    stat, p = mannwhitneyu(wine.iloc[:, i][wine['quality']<=5], wine.iloc[:, i][wine['quality']>5])
    print(names[i], ': p = %.3f' % (p))

# Se comprueba que para que cada prueba se obtiene un p-valor superior al nivel de signiﬁcación preﬁjado α = 0.05. Si esto se cumple, entonces se considera que dicha variable sigue una distribución normal.
#
# Por lo tanto aquellas variables con un p-valor mayor al nivel de significación son:
#
# **Residual sugar, pH**

# * #### **Test Kruskal-Wallis**

# En el caso que consideremos más de dos grupos, utilizamos el Test Kruskal-Wallis. Para ello utilizamos 3 niveles con la variable 'quality'

# In[27]:

from scipy.stats import kruskal

for i in range(len(names)):
    stat, p = kruskal(wine.iloc[:, i][wine['quality']<5],
                      wine.iloc[:, i][(wine['quality']==5) | (wine['quality']==6)],
                      wine.iloc[:, i][wine['quality']>6])
    print(names[i], ': p = %.3f' % (p))

# En este caso la única variable significativa, casi in extremis, es **residual sugar**

# ### ***Correlación***

# Procedemos a realizar un análisis de correlación entre las distintas variables para determinar cuáles de ellas ejercen una mayor inﬂuencia sobre la calidad del vino. Se utilizará el coeﬁciente de correlación de Person.

# In[28]:

# Matriz de correlaciones

f, ax = plt.subplots(figsize=(10, 6))
corr = wine.corr()
hm = sns.heatmap(round(corr,2), annot=True, ax=ax, cmap="Reds",fmt='.2f',
            linewidths=.05)
f.subplots_adjust(top=0.93)
t= f.suptitle('Wine Attributes Correlation Heatmap', fontsize=14)

# In[29]:

#correlations = wine.corr()['quality'].drop(['quality'])
corr['quality'].sort_values(ascending=False)

# Observamos que 'alcohol', 'sulphates', 'citric acid'  'fixed acidity' tienen una correlación máxima con 'quality'

# ### ***Regresiones***

# In[30]:

# definimos las variables X (atributos) y target('rating')
X = wine.drop(columns = ['quality','rating'])
y1 = wine['rating']

# In[31]:

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn import metrics

# Creamos los datos de entrenamiento y test utilizando train_test_split. Un 20% de los datos se utilizan para test y el 75%
# para entrenamiento
X_train, X_test, y_train, y_test = train_test_split(X, y1, test_size=0.2, random_state=1)

# In[32]:

# Ajustamos la regresión lineal a los datos de entrenamiento
regressor = LinearRegression()
regressor.fit(X_train,y_train)

# Con ello conseguimos los coeficientes de 10 atributos
regressor.coef_

# In[33]:

# Para predecir la calidad del vino con este módelo, usamos 'predict()'
train_pred = regressor.predict(X_train)
test_pred = regressor.predict(X_test)

# La precisión del Regresor Lineal

# In[34]:

accuracy = regressor.score(X_test, y_test)
print(accuracy)

# El RMSE es una medida que se utiliza para identificar las diferencias entre los valores.
# Si hemos construido un buen módelo el RMSE debe ser muy similar.
# Si el RMSE es mayor para el test que para en train, seguramente tengamos un mal ajuste de los datos.

# In[35]:

# Calculamos el  Root mean squared error (RMSE) para los datos de training como de test
train_rmse = mean_squared_error(train_pred, y_train) ** 0.5
print(train_rmse)
test_rmse = mean_squared_error(test_pred, y_test) ** 0.5
print(test_rmse)

# Los resultados para RMSE del train y test son bastante similares, por lo tanto hemos construido un buen módelo.

# In[36]:

print('Mean Absolute Error:', metrics.mean_absolute_error(y_test, test_pred))
print('Mean Squared Error:', metrics.mean_squared_error(y_test, test_pred))
print('RMSE:', np.sqrt(metrics.mean_squared_error(y_test, test_pred)))

# Mostramos los coeficientes para cada atributo
coeff = pd.DataFrame(regressor.coef_, list(data_wine.columns))
coeff.columns = ['Coefficient']
print(coeff)

# Estos coeficientes explican que si aumentamos en 1 unidad los 'sulphates', la calidad del vino aumentará un 0.50.
# Igualmente para otras variables. Las negativas como density 1 unidad de este atributo reducirá la calidad del vino en 15.

# ### ***Random Forest - Arbol de decisión***

# In[37]:

# llamamos al clasificador
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, log_loss
from sklearn.metrics import confusion_matrix

rfc = RandomForestClassifier(random_state=1)

# In[38]:

# Calculamos la  k-fold cross validation en el conjunto de entrenamiento para obtener la precisión

cv_scores = cross_val_score(rfc, X_train, y_train, cv=5, scoring='accuracy')

print('La precisión en cada iteración {}'.format(cv_scores))
print('La precisión media {}'.format(cv_scores.mean()))

# In[39]:

rfc.fit(X_train, y_train)
pred_rfc = rfc.predict(X_test)

# In[40]:

print(accuracy_score(y_test, pred_rfc))
print(log_loss(y_test, pred_rfc))

# In[41]:

print(confusion_matrix(y_test, pred_rfc))

# Hay un total de 61 errores de clasificación

# Una prediccion de 80% parece bastate buena para seleccionar el árbol de decisión como un clasificador recomendable para clasificar el vino tinto como '1' (>5) or "0" (<=5).

# # 5. Representación de los resultados a partir de tablas y gráficas.

# De los resultados anteriores, seleccionamos los atributos seleccionados en el apartado anterior:
#
# 'residual sugar', 'pH','volatile acidity','total sulfur dioxide', 'sulphates', 'alcohol','chlorides', 'citric acid', 'density'

# In[42]:

# Si consideramos algunas de las variables significativas encontradas en el apartado anterior

feats_to_explore = ['density','sulphates','volatile_acidity','citric_acid','chlorides', 'alcohol']

for i in feats_to_explore:
    plt.figure(figsize=(5,3))
    wine[i][wine['rating'] == '0'].plot.hist(alpha = 0.5, color = 'red', title = i)
    wine[i][wine['rating'] == '1'].plot.hist(alpha = 0.5, color = 'green')
    plt.show()

# ### ALCOHOL vs RATING

# In[43]:

bx = sns.boxplot(x="rating", y='alcohol', data = wine)
bx.set(xlabel='Rating Calidad', ylabel='Alcohol', title='Alcohol vs Rating Calidad')

# ### DENSITY vs RATING

# In[44]:

bx = sns.boxplot(x="rating", y='density', data = wine)
bx.set(xlabel='Rating', ylabel='Sulphates', title='Density vs Rating')

# ### SULPHATES vs RATING

# In[45]:

bx = sns.boxplot(x="rating", y='sulphates', data = wine)
bx.set(xlabel='Rating', ylabel='Sulphates', title='Sulphates vs Rating')

# ### CITRIC ACID vs RATING

# In[46]:

bx = sns.boxplot(x="rating", y='citric_acid', data = wine)
bx.set(xlabel='Rating', ylabel='Citric Acid', title='Citric Acid vs Rating')

# ### VOLATILE ACIDITY vs RATING

# In[47]:

x = sns.boxplot(x="rating", y='volatile_acidity', data = wine)
bx.set(xlabel='Rating', ylabel='Volatile Acidity', title='Volatile Acidity vs Rating')

# ### CHLORIDES vs RATING

# In[48]:

bx = sns.boxplot(x="rating", y="chlorides", data = wine);
bx.set(xlabel='Ratings', ylabel='chlorides', title='chlorides vs Ratings')

# ### LINEAR REGRESSION

# In[49]:

sns.lmplot(x = "alcohol", y = "residual_sugar", col = "rating", data = wine)

# In[50]:

sns.lmplot(x = "citric_acid", y = "volatile_acidity", col = "rating", data = wine)

# ### Biografia
#
# P. Cortez, A. Cerdeira, F. Almeida, T. Matos and J. Reis. Modeling wine preferences by data mining from physicochemical properties. In Decision Support Systems, Elsevier, 47(4):547-553, 2009
#
# https://github.com/ranjitkumarpatnaik/winequality-redwine/blob/master/winequality-red.ipynb
# https://www.kaggle.com/vishalyo990/prediction-of-quality-of-wine
# http://rpubs.com/jeswin_george/explore_red_quality_wines
# https://github.com/paulgx/tipologiapractica2/blob/master/practica2Tipologia.pdf
# https://www.kaggle.com/datacog314/tutorial-machine-learning-interpretability
#
# https://machinelearningmastery.com/a-gentle-introduction-to-normality-tests-in-python/
#
# Donde se pueden descargar los datos:
#
# https://archive.ics.uci.edu/ml/datasets/wine+quality

ext: py

reason: import numpy,from numpy,from scipy


text:
'''
Author: Mario Liu
Description: Module used to sample from depth matrices.
'''

import time
import numpy as np
from scipy import sparse

def createSamples(depth, perc_samples):
    '''
    Grabs samples from original data.
    Sparse identity matrices are used to conserve memory since identity
    matrices are mostly zeros.

    Args:
        depth: Depth matrix to grab samples from

        perc_samples: Percent sampling rate

    Returns:
        list: A list (of size less than or =
        perc_samples * # of cells in depth matrix)
        of flattened indices that correspond to non-NaN
        values in the given depth matrix

        list: The actual depth values at the indices
        given by the previous list
    '''
    '''
    Code adapted from sparse-depth-sensing/lib/sampling/createSamples.m
    '''
    if (perc_samples > 1.0) or (perc_samples < 0.0):
        print('perc_samples must be between 0 and 1')
        exit(1)

    height = depth.shape[0]
    width = depth.shape[1]
    N = height * width
    K = int(N * perc_samples)

    xGT = depth.flatten()
    # random sampling
    rand = np.random.permutation(N)[:K]
    # xGT[rand] permutes the xGT matrix according to the rand matrix, creates
    # a permuted matrix of size len(rand)
    a = np.isnan(xGT[rand])
    # rand[np.nonzero(~a)[0]] permutes the rand matrix according to the
    # nonzero(~a) matrix, creates a permuted matrix of size len(nonzero(~a))
    samples = rand[np.nonzero(~a)[0]]

    '''
    very slow way of getting depth values at indices in samples
    '''
    # Rfull = sparse.eye(N)
    # # for i in samples, takes the ith row from Rfull
    # R = Rfull.tocsr()[samples,:]
    # # scipy sparse matrix vector product, xGT is size n x 1
    # vec = R*xGT

    '''
    faster way of getting depth values at indices in samples
    '''
    vec = []
    for i in samples:
        vec.append(xGT[i])

    return samples, vec

def main():
    import matplotlib.pyplot as plt

    h = 12
    w = 16
    perc_samples = 0.5

    depth = np.zeros((h, w))
    depth.fill(np.nan)
    for _ in range(int((h * w) * 0.4)):
        y, x = int(h * np.random.sample()), int(w * np.random.sample())
        depth[y, x] = 6.0 * np.random.sample()
    nonNan = depth[~np.isnan(depth)]

    t1 = time.time()
    samples, vec = createSamples(depth, perc_samples)
    t2 = time.time()
    print('Sampling fraction: ' + str(perc_samples))
    print('Time to create samples: ' + str(t2 - t1))
    print('frac of samples: ' + str(len(samples)/float(len(nonNan))))
    print('samples: ' + str(samples))
    print('vec: ' + str(vec))

    plt.figure()

    plt.title('Samples ({0}% of All Pixels)'.format(perc_samples * 100))
    plt.imshow(depth, cmap='plasma')
    plt.colorbar()

    # mark each sample with red dot
    for samp in samples:
        x = samp % w
        y = samp / w
        plt.plot(x, y, 'wo', markersize=5)
        plt.plot(x, y, 'ro', markersize=3)

    plt.show()

if __name__ == "__main__":
    main()

max_stars_repo_name: marioliu/AutonomousQuadblade

max_issues_repo_name: marioliu/AutonomousQuadblade

max_forks_repo_name: marioliu/AutonomousQuadblade

reason: import numpy,from scipy


text:
import os
from skimage import io, transform
import torch
import torchvision
import argparse
from torch.autograd import Variable
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.data import Dataset, DataLoader
from torchvision import transforms
import numpy as np
from PIL import Image
import glob

from data_loader import RescaleT
from data_loader import ToTensorLab
from data_loader import SalObjDataset

from model import U2NET
from model import U2NETP

# 数据归一化
def normPRED(d):
    ma = torch.max(d)
    mi = torch.min(d)
    dn = (d-mi)/(ma-mi)

    return dn
# 参数
def get_args():
    parser = argparse.ArgumentParser(description='Train the U2net on images and target masks')
    parser.add_argument("--img_dir", help="data directory", default="/home/gaozijing/gzj/PRseg/data/dataset1/testuint8/", type=str)
    parser.add_argument("--modelname", help="model name", default="u2net", type=str)
    parser.add_argument('--scale', type=int, default=320, help='resize shape')
    parser.add_argument("--themodel", help="model name", default="u2net_bce_itr_106000_train_0.253333_tar_0.008061", type=str)
    return parser.parse_args()

def main():
    args = get_args()
    model_name= args.modelname
    image_dir = args.img_dir
    # 预测的保存路径
    prediction_dir = './test_data/' +args.themodel + '_results/'
    if not os.path.exists(prediction_dir):
        os.makedirs(prediction_dir)
    # 模型的路径
    model_dir = './saved_models/'+ model_name + '/' + args.themodel + '.pth'

    img_name_list = glob.glob(image_dir + '*')
    test_salobj_dataset = SalObjDataset(img_name_list = img_name_list,lbl_name_list = [],transform=transforms.Compose([RescaleT(args.scale),ToTensorLab(flag=0)]))
    test_salobj_dataloader = DataLoader(test_salobj_dataset,
                                        batch_size=1,
                                        shuffle=False,
                                        num_workers=1)

    net = U2NET(3,1)
    net.load_state_dict(torch.load(model_dir))
#     若多卡训练，则取消下一行的注释
#     net.load_state_dict({k.replace('module.',''):v for k,v in torch.load(model_dir).items()})
    if torch.cuda.is_available():
        net.cuda()
    net.eval()

    for i_test, data_test in enumerate(test_salobj_dataloader):
        print(i_test)
        inputs_test = data_test['image']
        inputs_test = inputs_test.type(torch.FloatTensor)

        if torch.cuda.is_available():
            inputs_test = Variable(inputs_test.cuda())
        else:
            inputs_test = Variable(inputs_test)

        d1,d2,d3,d4,d5,d6,d7= net(inputs_test)
        # 顶层输出
        pred = d1[:,0,:,:]
        # 输出归一化
        pred = normPRED(pred)
        del d1,d2,d3,d4,d5,d6,d7
        image_name = img_name_list[i_test]
        predict = pred
        predict = predict.squeeze()
        predict_np = predict.cpu().data.numpy()
        # 转换为uint8的彩图格式
        im = Image.fromarray(predict_np*255).convert('RGB')
        img_name = image_name.split("/")[-1]
        # 放缩至原图大小
        imo = im.resize((628,628),resample=Image.BILINEAR)
        pb_np = np.array(imo)
        # 获取文件名
        aaa = img_name.split(".")
        bbb = aaa[0:-1]
        imidx = bbb[0]
        for i in range(1,len(bbb)):
            imidx = imidx + "." + bbb[i]
        imo.save(prediction_dir+imidx+'.png')

if __name__ == "__main__":
    main()

max_stars_repo_name: ZjGaothu/Cell-Instance-Segmentation

max_issues_repo_name: ZjGaothu/Cell-Instance-Segmentation

max_forks_repo_name: ZjGaothu/Cell-Instance-Segmentation

reason: import numpy


text:
from os import path
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import os

from wordcloud import WordCloud, STOPWORDS, ImageColorGenerator

# get data directory (using getcwd() is needed to support running example in generated IPython notebook)
d = path.dirname(__file__) if "__file__" in locals() else os.getcwd()

# Read the whole text.
text = open(path.join(d, 'dialogues.txt')).read()

mask = np.array(Image.open(path.join(d, "images/steve_jobs.jpeg")))
stopwords = set(STOPWORDS)

wc = WordCloud(width=1000, height=600, background_color="black", max_words=2000, mask=mask,
               stopwords=stopwords, max_font_size=90, random_state=42)
# generate word cloud
wc.generate(text)

# create coloring from image
image_colors = ImageColorGenerator(mask)

# show
# fig, axes = plt.subplots(1, 3)
# axes[0].imshow(wc, interpolation="bilinear")
# # recolor wordcloud and show
# # we could also give color_func=image_colors directly in the constructor
# axes[1].imshow(wc.recolor(color_func=image_colors), interpolation="bilinear")
# axes[2].imshow(alice_coloring, cmap=plt.cm.gray, interpolation="bilinear")
# for ax in axes:
#     ax.set_axis_off()
# plt.show()

plt.figure(figsize=[7,7])
plt.imshow(wc.recolor(color_func=image_colors), interpolation="bilinear")
plt.axis("off")
_=plt.show()

max_stars_repo_name: FOSSNITRR/word-cloud-generator

max_issues_repo_name: FOSSNITRR/word-cloud-generator

max_forks_repo_name: FOSSNITRR/word-cloud-generator


text:
import numpy as np
from matplotlib import pyplot as plt
from scipy import sparse

eye = np.eye(4)
print("NumPy array:\n{}".format(eye))

sparse_matrix = sparse.csr_matrix(eye)
print("\nSciPy sparse CSR matrix:\n{}".format(sparse_matrix))

x = np.linspace(-10, 10, 100)
y = np.sin(x)
plt.plot(x, y)
plt.show()


text:
import sys
import cv2
import numpy as np
import time
import imutils
from matplotlib import pyplot as plt

# Functions
import HSV_filter as hsv
import shape_recognition as shape
import triangulation as tri
import no_mask
import predict_real_time as pre
#import calibration as calib

# Open both cameras
frame_left = cv2.imread('/home/daniel/Escritorio/MSc AAI/IRP/blender projects/RANGER/Render/cube/Stereo Vision/cube_stereo_006_left.png')
frame_right = cv2.imread('/home/daniel/Escritorio/MSc AAI/IRP/blender projects/RANGER/Render/cube/Stereo Vision/cube_stereo_006_right.png')

B = 10.9878               #Distance between the cameras [m]
f = 25               #Camera lense's focal length [mm]
alpha = 71.5078        #Camera field of view in the horizontal plane [degrees]

################## CALIBRATION #########################################################

#frame_right, frame_left = calib.undistorted(frame_right, frame_left)

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

# # APPLYING HSV-FILTER:
# mask_right = hsv.add_HSV_filter(frame_right)
# mask_left = hsv.add_HSV_filter(frame_left)

# # Result-frames after applying HSV-filter mask
# res_right = cv2.bitwise_and(frame_right, frame_right, mask=mask_right)
# res_left = cv2.bitwise_and(frame_left, frame_left, mask=mask_left)

# # APPLYING SHAPE RECOGNITION:
# circles_right = shape.find_circles(frame_right, mask_right)
# circles_left  = shape.find_circles(frame_left, mask_left)

# APPLYING SHAPE RECOGNITION WITHOUT THE MASK:
circles_right = no_mask.find_circles(frame_right)
circles_left  = no_mask.find_circles(frame_left)

# Hough Transforms can be used aswell or some neural network to do object detection
# bounding_box_right = pre.find_bounding_box(frame_right)
# bounding_box_left  = pre.find_bounding_box(frame_left)

        ################## CALCULATING BALL DEPTH #########################################################

# If no ball can be caught in one camera show text "TRACKING LOST"
if np.all(circles_right) == None or np.all(circles_left) == None:
    cv2.putText(frame_right, "TRACKING LOST", (75,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,0,255),2)
    cv2.putText(frame_left, "TRACKING LOST", (75,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,0,255),2)

else:
# Function to calculate depth of object. Outputs vector of all depths in case of several balls.
# All formulas used to find depth is in video presentaion
    depth = tri.find_depth(circles_right, circles_left, frame_right, frame_left, B, f, alpha)

    # cv2.putText(frame_right, "TRACKING", (75,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (124,252,0),2)
    # cv2.putText(frame_left, "TRACKING", (75,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (124,252,0),2)
    # cv2.putText(frame_right, "Distance: " + str(round(depth,3)), (200,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (124,252,0),2)
    # cv2.putText(frame_left, "Distance: " + str(round(depth,3)), (200,50), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (124,252,0),2)
    # Multiply computer value with 205.8 to get real-life depth in [cm]. The factor was found manually.
    cv2.putText(frame_right, "{0:.2f}m".format(depth),
        (frame_right.shape[1] - 400, frame_right.shape[0] - 20), cv2.FONT_HERSHEY_SIMPLEX,
        2.0, (0, 255, 0), 3)
    cv2.putText(frame_left, "{0:.2f}m".format(depth),
        (frame_left.shape[1] - 400, frame_left.shape[0] - 20), cv2.FONT_HERSHEY_SIMPLEX,
        2.0, (0, 255, 0), 3)
    print("Depth: ", depth)

# Resize the frames
frame_right = cv2.resize(frame_right, (960, 540))
frame_left = cv2.resize(frame_left, (960, 540))
# mask_right = cv2.resize(mask_right, (960, 540))
# mask_left = cv2.resize(mask_left, (960, 540))

# Show the frames
# cv2.imshow("frame right", frame_right)
# cv2.imshow("frame left", frame_left)
# # cv2.imshow("mask right", mask_right)
# # cv2.imshow("mask left", mask_left)
# cv2.waitKey(0)

# cv2.destroyAllWindows()

cv2.imwrite('/home/daniel/Escritorio/MSc AAI/IRP/Results/Stereo Vision/cube_stereo_006_left.png', frame_left)
cv2.imwrite('/home/daniel/Escritorio/MSc AAI/IRP/Results/Stereo Vision/cube_stereo_006_right.png', frame_right)


text:
# -*- coding: utf-8 -*-
import os
import imgaug as ia
from imgaug import augmenters as iaa
import scipy
import shutil

path = 'C:\\Users\\HEVEN\\Desktop\\background_ilseok\\background'

augmentation_path = 'C:\\Users\\HEVEN\\Desktop\\background_ilseok\\augmentation'

for root, dirs, files in os.walk(path):
    for file in files:
        if file.endswith(".jpg"):
            for light in range(10):
                jpg_path = path + "\\" + file
                txt_path = path + "\\" + file[:-4] + ".txt"

                image = scipy.misc.imread(jpg_path, flatten=False, mode="RGB")

                if light < 5:
                    seq = iaa.Sequential([
                            iaa.Multiply( 1 - 0.15 * (light) )
                    ])
                else:
                    seq = iaa.Sequential([
                            iaa.Multiply( 1 + 0.3 * (light-4) )
                    ])

                image_aug = seq.augment_images([image])[0]

                new_jpg_path = augmentation_path + "\\new_" + file[:-4] + "_%d"%light + ".jpg"
                new_txt_path = augmentation_path + "\\new_" + file[:-4] + "_%d"%light + ".txt"

                scipy.misc.imsave(new_jpg_path, image_aug)
                shutil.copy2(txt_path, new_txt_path)

max_stars_repo_name: HongBeenKim/pams-skku

max_issues_repo_name: HongBeenKim/pams-skku

max_forks_repo_name: HongBeenKim/pams-skku

reason: import scipy


text:
import numpy as np
# 把数据特征标准化
def feature_nomalize(X):
	mu = np.mean(X,0)
	sigma = np.std(X,0,ddof=1)
	X_norm = (X - mu)/sigma

	return X_norm,mu,sigma


text:
#!/usr/bin/env python
# By Myles McKay
# June 7, 2020
# ZTF crossmatch with X-Ray Binaries (ROSAT catalog)

import io
import time
import datetime
import argparse
import logging
import gzip
from threading import Lock
from copy import deepcopy
import numpy as np
import pandas as pd
import fastavro
from astropy.coordinates import SkyCoord, match_coordinates_sky
import astropy.units as u
from astropy.io import fits
from aiokafka import AIOKafkaConsumer
import asyncio
from concurrent.futures import ThreadPoolExecutor
import functools

from astroquery.simbad import Simbad
from .constants import LOGGING, BASE_DIR, DB_DIR, SIMBAD_EXCLUDES # FITS_DIR
from .db_caching import create_connection, insert_data, update_value, insert_lc_dataframe, get_cached_ids

# Example command line execution:

# for date 210119, program_id 1, 64 cores, suffix 64t
# $ python ztf_rosat_crossmatch.py 210119 1 64 64t

SIGMA_TO_95pctCL = 1.95996
FITS_DIR = "../local/cutouts_debug"

def exception_handler(func):
    def wrapper(*args, **kwargs):
        try:
            return func(*args, **kwargs)
        except Exception as e:
            logging.exception(e)

    return wrapper

@exception_handler
def read_avro_file(fname):
    """Reads a single packet from an avro file stored with schema on disk."""
    with open(fname, "rb") as f:
        freader = fastavro.reader(f)
        for packet in freader:
            return packet

@exception_handler
def read_avro_bytes(buf):
    """Reads a single packet from an avro file stored with schema on disk."""
    with io.BytesIO(buf) as f:
        freader = fastavro.reader(f)
        for packet in freader:
            return packet

@exception_handler
def get_candidate_info(packet):
    return {"ra": packet["candidate"]["ra"], "dec": packet["candidate"]["dec"],
            "object_id": packet["objectId"], "candid": packet["candid"]}

@exception_handler
def save_cutout_fits(packet, output):
    """Save fits cutouts from packed into output."""
    objectId = packet["objectId"]
    pid = packet["candidate"]["pid"]
    for im_type in ["Science", "Template", "Difference"]:
        with gzip.open(io.BytesIO(packet[f"cutout{im_type}"]["stampData"]), "rb") as f:
            with fits.open(io.BytesIO(f.read())) as hdul:
                hdul.writeto(f"{output}/{objectId}_{pid}_{im_type}.fits", overwrite=True)

@exception_handler
def make_dataframe(packet, repeat_obs=True):
    """Extract relevant lightcurve data from packet into pandas DataFrame."""
    df = pd.DataFrame(packet["candidate"], index=[0])
    if repeat_obs:
        df["ZTF_object_id"] = packet["objectId"]
        return df[["ZTF_object_id", "jd", "fid", "magpsf", "sigmapsf", "diffmaglim"]]

    df_prv = pd.DataFrame(packet["prv_candidates"])
    df_merged = pd.concat([df, df_prv], ignore_index=True)
    df_merged["ZTF_object_id"] = packet["objectId"]
    return df_merged[["ZTF_object_id", "jd", "fid", "magpsf", "sigmapsf", "diffmaglim"]]

@exception_handler
def load_rosat():
    # Open ROSAT catalog
    rosat_fits = fits.open("/epyc/data/rosat_2rxs/cat2rxs.fits")

    # Make ROSAT data into a pandas dataframe
    rosat_data = rosat_fits[1].data
    dfx = pd.DataFrame(rosat_data)
    rosat_fits.close()

    # exclude sources that are not observable
    dfx = dfx[dfx.DEC_DEG >= -30]

    # List ROSAT error position [arcsec]
    dfx["err_pos_arcsec"] = np.sqrt(((dfx.XERR * 45) ** 2. + (dfx.YERR * 45) ** 2.) + 0.6 ** 2.) * SIGMA_TO_95pctCL

    # Put ROSAT ra and dec list in SkyCoord [degrees]
    rosat_skycoord = SkyCoord(ra=dfx.RA_DEG, dec=dfx.DEC_DEG,
                              frame="icrs", unit=(u.deg))

    return dfx[["IAU_NAME", "RA_DEG", "DEC_DEG", "err_pos_arcsec"]], rosat_skycoord

@exception_handler
def not_moving_object(packet):
    """Check if there are > 2 detections separated by at least 30 minutes.

    Parameters:
       avro_packet: dictionary
          Extracted data from the avro file

    Return:
        real: boolean
            True if there is another detection > 30 minutes from the triggering detection
    """

    date0 = float(packet['candidate']['jd'])
    if packet['prv_candidates'] is None:
        logging.info("prv_candidates is None")
        return False
    if len(packet['prv_candidates']) == 0:
        logging.debug("No previous detections of {packet['objectId']}")
        return False

    for prv_candidate in packet['prv_candidates']:
        if prv_candidate['candid'] is not None:
            date1 = float(prv_candidate['jd'])
            diff_date = date0 - date1
            if diff_date < (0.5 / 24):
                continue
            else:
                # logging.debug(f"Previous detection of {packet['objectId']} {diff_date} days earlier")
                return True

    return False

customSimbad = Simbad()
customSimbad.add_votable_fields("otype(3)")
customSimbad.add_votable_fields("otypes(3)")  # returns a '|' separated list of all the otypes

@exception_handler
def query_simbad(ra, dec):
    sc = SkyCoord(ra, dec, unit=u.degree)
    result_table = customSimbad.query_region(sc, radius=2 * u.arcsecond)
    return result_table

def is_excluded_simbad_class(ztf_source):
    """Is the object in Simbad, with object types we reject (e.g., AGN)?
    """
    try:
        result_table = query_simbad(ztf_source['ra'], ztf_source['dec'])
        if result_table is None:
            logging.info(f"{ztf_source['object_id']} not found in Simbad")
            return False  # not in Simbad
        else:
            # for now let's just take the closest match
            otype = result_table['OTYPE_3'][0].decode("utf-8")
            simbad_id = result_table['MAIN_ID'][0].decode("utf-8")
            if otype in SIMBAD_EXCLUDES:
                logging.info(f"{ztf_source['object_id']} found in Simbad as {simbad_id} ({otype}); omitting")
                return True
            else:
                logging.info(f"{ztf_source['object_id']} found in Simbad as {simbad_id} ({otype}); saving")
                return False

    except Exception as e:
        # if this doesn't work, record the exception and continue
        logging.exception(f"Error querying Simbad for {ztf_source['object_id']}", e)
        return False

def ztf_rosat_crossmatch(ztf_source, rosat_skycoord, dfx):
    """
    Cross match ZTF and ROSAT data using astropy.coordinates.SkyCoord

    Parameters:
                ztf_source: dict
                    {'ra': float (degrees), 'dec': float (degrees),
                    'object_id': string, 'candid': int}

                rosat_skycoord: astropy.coordinates.SkyCoord
                    ROSAT catalog in astropy.coordinates.SkyCoord

                dfx: pandas.DataFrame
                    slim ROSAT catalog

    Return:
            matched_source: dict or None
                if a source matches, return
                    {'ra': float (degrees), 'dec': float (degrees),
                    'source_name': string, 'sep2d': float (arcsec)}
                else None
    """
    try:
        # Input avro data ra and dec in SkyCoords
        avro_skycoord = SkyCoord(ra=ztf_source["ra"], dec=ztf_source["dec"],
                                 frame="icrs", unit=(u.deg))

        # Finds the nearest ROSAT source's coordinates to the avro files ra[deg] and dec[deg]
        match_idx, match_sep2d, _ = avro_skycoord.match_to_catalog_sky(rosat_skycoord)

        match_row = dfx.iloc[match_idx]

        matched = match_sep2d[0] <= match_row["err_pos_arcsec"] * u.arcsecond

        match_result = {"match_name": match_row["IAU_NAME"],
                        "match_ra": match_row["RA_DEG"],
                        "match_dec": match_row["DEC_DEG"],
                        "match_err_pos": match_row["err_pos_arcsec"],
                        "match_sep": match_sep2d[0].to(u.arcsecond).value}

        if matched:
            logging.info(
                f"{ztf_source['object_id']} ({avro_skycoord.to_string('hmsdms')}; {ztf_source['candid']}) matched {match_result['match_name']} ({match_result['match_sep']:.2f} arcsec away)")
            return match_result

        else:
            # logging.debug(f"{ztf_source['object_id']} ({avro_skycoord.to_string('hmsdms')}) did not match (nearest source {match_result['match_name']}, {match_result['match_sep']:.2f} arcsec away")
            return None
    except Exception as e:
        logging.exception(f"Unable to crossmatch {ztf_source['object_id']} with ROSAT", e)

@exception_handler
def save_to_db(packet, otype, sources_saved, lock_sources_saved, database, interest):
    """Save matches to database
    """
    ztf_object_id = packet["objectId"]
    data_to_update = {"SIMBAD_otype": f'"{otype}"', "ra": packet["candidate"]["ra"],
                      "dec": packet["candidate"]["dec"], "SIMBAD_include": interest}
    logging.info(f"Saving new source {ztf_object_id} to database.")

    with lock_sources_saved:
        conn = create_connection(database)
        if ztf_object_id in sources_saved:
            logging.info(f"{ztf_object_id} already saved in time between simbad check and now")
            dflc = make_dataframe(packet, repeat_obs=True)
        else:
            update_value(conn, data_to_update, f'ZTF_object_id = "{ztf_object_id}"')
            dflc = make_dataframe(packet, repeat_obs=False)
            # cache_ZTF_object(conn, (ztf_object_id, simbad_otype, ra, dec, rosat_iau_name))
            logging.debug(f"Successfully saved new source {ztf_object_id} to database.")
        insert_lc_dataframe(conn, dflc)
        logging.debug(f"Successfully ingested lightcurve data from {ztf_object_id} to database.")
        save_cutout_fits(packet, FITS_DIR)
        logging.debug(f"Successfully saved cutouts of {ztf_object_id}")
        sources_saved.update((ztf_object_id,))
        conn.close()

@exception_handler
def check_for_new_sources(packets_to_simbad, sources_saved, lock_sources_saved, database):
    """Checks the packets_to_simbad for ZTF objects not previously saved to the database.
    """
    logging.debug(f"Checking for new sources...")
    with lock_sources_saved:
        new_packets = [packet for packet in packets_to_simbad if packet["objectId"] not in sources_saved]
        old_packets = [packet for packet in packets_to_simbad if packet["objectId"] in sources_saved]
        logging.debug("New sources: {}".format(", ".join([packet["objectId"] for packet in new_packets])))
        if len(new_packets) < len(packets_to_simbad):
            logging.info(f"{len(packets_to_simbad) - len(new_packets)} seen before")

    logging.debug("Updating info for previously seen sources...")
    for packet in old_packets:
        ztf_object_id = packet["objectId"]
        conn = create_connection(database)
        dflc = make_dataframe(packet, repeat_obs=True)
        insert_lc_dataframe(conn, dflc)
        logging.debug(f"Successfully updated lightcurve data from {ztf_object_id} to database.")
        save_cutout_fits(packet, FITS_DIR)
        logging.debug(f"Successfully updated cutouts of {ztf_object_id}")
        conn.close()

    return new_packets

@exception_handler
def process_packet(packet, rosat_skycoord, dfx, saved_packets, lock, sources_seen, lock_sources_seen, database):
    """Examine packet for matches in the ROSAT database. Save object to database if match found"""
    if packet["candidate"]["drb"] < 0.8:  # if packet real/bogus score is low, ignore
        return
    ztf_source = get_candidate_info(packet)
    # conn = create_connection(database)
    with lock:
        with lock_sources_seen:
            if packet["objectId"] in sources_seen:
                # logging.debug(f"{packet['objectId']} already known match, adding packet to packets_from_kafka")
                saved_packets.append(packet)
                sources_seen.update((packet["objectId"],))
                # conn.close()
                logging.debug(f"Total of {len(saved_packets)} saved for query.")
            else:
                matched_source = ztf_rosat_crossmatch(ztf_source, rosat_skycoord, dfx)
                if (matched_source is not None) and not_moving_object(packet):
                    # logging.debug("adding packet to packets_from_kafka")
                    saved_packets.append(packet)
                    sources_seen.update((packet["objectId"],))
                    conn = create_connection(database)
                    data_to_insert = {"ZTF_object_id": packet["objectId"],
                                      "ROSAT_IAU_NAME": matched_source["match_name"]}
                    insert_data(conn, "ZTF_objects", data_to_insert)
                    logging.debug(f"Successfully saved {packet['objectId']} to database")
                    conn.close()
                    logging.debug(f"Total of {len(saved_packets)} saved for query.")

@exception_handler
def check_simbad_and_save(packets_to_simbad, lock_packets_to_simbad, sources_saved, lock_sources_saved, database):
    with lock_packets_to_simbad:
        logging.info("Checking packets for new sources")
        new_packets_to_simbad = check_for_new_sources(packets_to_simbad, sources_saved,
                                                      lock_sources_saved, database)
        logging.debug(f"{len(packets_to_simbad) - len(new_packets_to_simbad)} sources already cached.")

        # Return if sources were previously seen and recorded
        if len(new_packets_to_simbad) == 0:
            logging.debug(f"All {len(packets_to_simbad)} were previously seen.")
            return None

        logging.info(f"{len(new_packets_to_simbad)} being sent to Simbad for matching")
        ras = []
        decs = []
        for packet in new_packets_to_simbad:
            ras.append(packet["candidate"]["ra"])
            decs.append(packet["candidate"]["dec"])

        sc = SkyCoord(ras, decs, unit=u.degree)

        MATCH_RADIUS = 2 * u.arcsecond
        try:
            result_table = customSimbad.query_region(sc, radius=MATCH_RADIUS)  # TODO: examine result_table for errors
            # annoyingly, the result table is unsorted and unlabeled, so we have
            # to re-match it to our packets
            if result_table is not None:
                sc_simbad = SkyCoord(result_table["RA"], result_table["DEC"],
                                     unit=("hourangle", "degree"))
                idx, sep2d, _ = match_coordinates_sky(sc, sc_simbad)
                # assert(len(sc) == len(idx))
            else:
                logging.info("No matches found in SIMBAD")
                return

        except Exception as e:
            logging.exception("Error querying Simbad", e)

        for i, packet in enumerate(new_packets_to_simbad):
            # check if we have a simbad match for each packet we sent
            if sep2d[i] <= MATCH_RADIUS:
                matched_row = result_table[idx[i]]
                otype = matched_row["OTYPE_3"].decode("utf-8")
                simbad_id = matched_row["MAIN_ID"].decode("utf-8")
                if otype in SIMBAD_EXCLUDES:
                    logging.info(f"{packet['objectId']} found in Simbad as {simbad_id} ({otype}); omitting")
                    save_to_db(packet, otype, sources_saved, lock_sources_saved, database, interest=0)
                else:
                    logging.info(f"{packet['objectId']} found in Simbad as {simbad_id} ({otype}); saving")
                    save_to_db(packet, otype, sources_saved, lock_sources_saved, database, interest=1)

            else:
                # no match in simbad,
                logging.info(f"{packet['objectId']} not found in Simbad")
                save_to_db(packet, None, sources_saved, lock_sources_saved, database, interest=0)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("date", type=str, help="UTC date as YYMMDD")
    parser.add_argument("program_id", type=int, help="Program ID (1 or 2)")
    parser.add_argument("threads", type=int, help="Number of threads (1 to 64)")
    parser.add_argument("suffix", type=str, help="suffix of db")

    args = parser.parse_args()

    if len(args.date) != 6:
        raise ValueError(f"Date must be specified as YYMMDD.  Provided {args.date}")

    if args.program_id not in [1, 2]:
        raise ValueError(f"Program id must be 1 or 2.  Provided {args.program_id}")

    if (args.threads < 1) or (args.threads > 64):
        raise ValueError(f"threads must be between 1 and 64.  Provided {args.threads}")

    database = DB_DIR + f"sqlite{args.suffix}.db"
    logging.debug(f"Database at {database}")

    kafka_topic = f"ztf_20{args.date}_programid{args.program_id}"
    kafka_server = "partnership.alerts.ztf.uw.edu:9092"

    now = datetime.datetime.now().strftime("%d%m%y_%H%M%S")
    LOGGING["handlers"]["logfile"]["filename"] = f"{BASE_DIR}/../logs/{kafka_topic}_{now}.log"
    logging.config.dictConfig(LOGGING)

    logging.debug(f"Args parsed and validated: {args.date}, {args.program_id}")

    # load X-ray catalogs
    dfx, rosat_skycoord = load_rosat()

    logging.info(f"Connecting to Kafka topic {kafka_topic}")

    loop = asyncio.get_event_loop()

    async def consume():
        consumer = AIOKafkaConsumer(
            kafka_topic,  # other_topic,
            loop=loop, bootstrap_servers=kafka_server,
            auto_offset_reset="earliest",
            value_deserializer=read_avro_bytes,
            group_id=f"uw_xray_{args.suffix}")
        # Get cluster layout and join group `my-group`
        await consumer.start()
        tstart = time.perf_counter()
        tbatch = tstart
        i = 0
        nmod = 1000
        pool = ThreadPoolExecutor(max_workers=args.threads)
        packets_from_kafka = []
        packets_to_simbad = []
        try:
            conn = create_connection(database)
        except Exception as e:
            logging.exception(f"Could not connect to {database}", e)
            return
        sources_saved = set(get_cached_ids(conn).values)
        sources_seen = set(get_cached_ids(conn).values)
        logging.info(f"{len(sources_saved)} sources previously seen")
        n0_sources = len(sources_saved)
        conn.close()
        lock_packets_from_kafka = Lock()
        lock_packets_to_simbad = Lock()
        lock_sources_saved = Lock()
        lock_sources_seen = Lock()
        logging.debug("begin ingesting messages")
        try:
            # Consume messages
            async for msg in consumer:
                i += 1
                if i % nmod == 0:
                    elapsed = time.perf_counter() - tstart
                    logging.info(f"Consumed {i} messages in {elapsed:.1f} sec ({i / elapsed:.1f} messages/s)")
                    with lock_sources_saved:
                        logging.info(f"Matched {len(sources_saved) - n0_sources} sources seen in {elapsed:.1f} sec")

                # query simbad in batches
                if time.perf_counter() - tbatch >= 40:
                    with lock_packets_from_kafka:
                        with lock_packets_to_simbad:
                            logging.debug("copying {} packets from kafka to simbad".format(len(packets_from_kafka)))
                            packets_to_simbad = deepcopy(packets_from_kafka)
                            packets_from_kafka = []
                    if len(packets_to_simbad) > 0:
                        logging.debug(f"{len(packets_to_simbad)} packets to query")
                        loop.run_in_executor(pool,
                                             functools.partial(check_simbad_and_save,
                                                               packets_to_simbad,
                                                               lock_packets_to_simbad,
                                                               sources_saved,
                                                               lock_sources_saved,
                                                               database))
                    tbatch = time.perf_counter()

                # print("consumed: ", msg.topic, msg.partition, msg.offset,
                #      msg.key, msg.timestamp)

                packet = msg.value

                loop.run_in_executor(pool,
                                     functools.partial(process_packet,
                                                       packet, rosat_skycoord, dfx, packets_from_kafka,
                                                       lock_packets_from_kafka,
                                                       sources_seen, lock_sources_seen, database))

        finally:
            # Will leave consumer group; perform autocommit if enabled.
            await consumer.stop()
            pool.shutdown(wait=True)
            # TODO: flush out saved packets

    loop.run_until_complete(consume())

reason: import numpy,import astropy,from astropy


text:
import numpy as np
import os
from mpEntropy import mpSystem
import matplotlib as mpl
from matplotlib.pyplot import cm
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from scipy.stats import norm
from scipy.signal import savgol_filter
from scipy.optimize import curve_fit
# This is a workaround until scipy fixes the issue
import warnings

warnings.filterwarnings(action="ignore", module="scipy", message="^internal gelsd")

# Load sysVar
sysVar = mpSystem("../interact_0.ini", plotOnly=True)

# Create plot folder
pltfolder = "./epsplots/"
if not os.path.exists(pltfolder):
    os.mkdir(pltfolder)

print("Plotting", end='')
mpl.use('Agg')
# minimum and maximum times to plot
min_time = 0
max_time = 3
inlay_min_time = 10
inlay_max_time = 100
inlay_log_min_time = 0
inlay_log_max_time = 3
# styles and stuff
avgstyle = 'dashed'
avgsize = 0.6
expectstyle = 'solid'
expectsize = 1
legend_size = 10
font_size = 10
# https://scipy.github.io/old-wiki/pages/Cookbook/Matplotlib/LaTeX_Examples.html
fig_width_pt = 246.0  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0 / 72.27  # Convert pt to inches
golden_mean = (np.sqrt(5) - 1.0) / 2.0  # Aesthetic ratio
fig_width = fig_width_pt * inches_per_pt  # width in inches
fig_height = fig_width * golden_mean  # height in inches
fig_size = [fig_width, fig_height]
# padding in units of fontsize
padding = 0.32

params = {
    'axes.labelsize': 10,
    'font.size': 10,
    'legend.fontsize': 10,
    'xtick.labelsize': 8,
    'ytick.labelsize': 8,
    'lines.linewidth': 1,
    'figure.figsize': fig_size,
    'legend.frameon': False,
    'legend.loc': 'best',
    'mathtext.default': 'rm'  # see http://matplotlib.org/users/customizing.html
}
plt.rcParams['agg.path.chunksize'] = 0
plt.rcParams.update(params)
plt.rc('text', usetex=True)
plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Arial']})

loavgpercent = sysVar.plotLoAvgPerc  # percentage of time evolution to start averaging
loavgind = int(loavgpercent * sysVar.dataPoints)  # index to start at when calculating average and stddev
loavgtime = np.round(loavgpercent * (sysVar.deltaT * sysVar.steps * sysVar.plotTimeScale), 2)

# stuff for averaging
if sysVar.boolPlotAverages:
    print(' with averaging from Jt=%.2f' % loavgtime, end='')
fwidth = sysVar.plotSavgolFrame
ford = sysVar.plotSavgolOrder

ent_array = np.loadtxt('../data/entropy.txt')

# multiply step array with time scale
step_array = ent_array[:, 0] * sysVar.plotTimeScale
min_index = int(min_time / step_array[-1] * len(step_array))
max_index = int(max_time / step_array[-1] * len(step_array))
inlay_min_index = int(inlay_min_time / step_array[-1] * len(step_array))
inlay_max_index = int(inlay_max_time / step_array[-1] * len(step_array))
inlay_log_min_index = int(inlay_log_min_time / step_array[-1] * len(step_array))
inlay_log_max_index = int(inlay_log_max_time / step_array[-1] * len(step_array))
#### Complete system Entropy
if os.path.isfile('../data/total_entropy.txt'):
    totent_array = np.loadtxt('../data/total_entropy.txt')
    plt.plot(totent_array[min_index:max_index, 0] * sysVar.plotTimeScale, totent_array[min_index:max_index, 1] * 1e13,
             linewidth=0.6, color='r')

    plt.grid()
    plt.xlabel(r'$J\,t$')
    plt.ylabel(r'Total system entropy $/ 10^{-13}$')
    plt.tight_layout(padding)
    ###
    plt.savefig(pltfolder + 'entropy_total.eps', format='eps', dpi=1000)
    plt.clf()
    print('.', end='', flush=True)

### Subsystem Entropy
fldat = open(pltfolder + 'ent_fluctuation_N' + str(sysVar.N) + '.txt', 'w')
fldat.write('N_tot: %i\n' % sysVar.N)
avg = np.mean(ent_array[loavgind:, 1], dtype=np.float64)
stddev = np.std(ent_array[loavgind:, 1], dtype=np.float64)
fldat.write('ssent_average: %.16e\n' % avg)
fldat.write('ssent_stddev: %.16e\n' % stddev)
fldat.write('ssent_rel._fluctuation: %.16e\n' % (stddev / avg))
fldat.close()

plt.plot(step_array[min_index:max_index], ent_array[min_index:max_index, 1], color='r')
if sysVar.boolPlotAverages:
    tavg = savgol_filter(ent_array[:, 1], fwidth, ford)
    plt.plot(step_array[loavgind:], tavg[loavgind:], linewidth=avgsize, linestyle=avgstyle, color='black')
plt.xlabel(r'$J\,t$')
plt.ylabel(r'$S\textsubscript{sys}$')
plt.tight_layout(padding)
plt.savefig(pltfolder + 'entropy_subsystem.eps', format='eps', dpi=1000)
plt.clf()

# Subsystem entropy with logarithmic inlay
plt.plot(step_array[min_index:max_index], ent_array[min_index:max_index, 1], color='r')
if sysVar.boolPlotAverages:
    tavg = savgol_filter(ent_array[:, 1], fwidth, ford)
    plt.plot(step_array[loavgind:], tavg[loavgind:], linewidth=avgsize, linestyle=avgstyle, color='black')
plt.xlabel(r'$J\,t$')
plt.ylabel(r'$S\textsubscript{sys}$')
a = plt.axes([.5, .35, .4, .4])
plt.semilogy(step_array[inlay_log_min_index:inlay_log_max_index],
             np.abs(avg - ent_array[inlay_log_min_index:inlay_log_max_index, 1]), color='r')
plt.ylabel(r'$|\,\overline{S}\textsubscript{sys} - S\textsubscript{sys}(t)|$')
plt.yticks([])
plt.savefig(pltfolder + 'entropy_subsystem_inlay_log.eps', format='eps', dpi=1000)
plt.clf()

# Subsystem entropy with inlay
plt.plot(step_array[min_index:max_index], ent_array[min_index:max_index, 1], color='r')
if sysVar.boolPlotAverages:
    tavg = savgol_filter(ent_array[:, 1], fwidth, ford)
    plt.plot(step_array[loavgind:], tavg[loavgind:], linewidth=avgsize, linestyle=avgstyle, color='black')
plt.xlabel(r'$J\,t$')
plt.ylabel(r'$S\textsubscript{sys}$')
a = plt.axes([.45, .35, .4, .4])
plt.plot(step_array[inlay_min_index:inlay_max_index], avg - ent_array[inlay_min_index:inlay_max_index, 1],
         linewidth=0.2, color='r')
plt.ylabel(r'$\overline{S}\textsubscript{sys} - S\textsubscript{sys}(t)$')
a.yaxis.tick_right()
tmp_ticks = list(a.get_xticks())
tmp_ticks.pop(0)
if tmp_ticks[-1] >= inlay_max_time:
    tmp_ticks.pop(-1)
a.set_xticks(tmp_ticks + [inlay_min_time])
plt.savefig(pltfolder + 'entropy_subsystem_inlay.eps', format='eps', dpi=1000)
plt.clf()
print('.', end='', flush=True)

# histogram of fluctuations
n, bins, patches = plt.hist(ent_array[loavgind:, 1] - avg, 51, normed=1, rwidth=0.8, align='mid')
(mu, sigma) = norm.fit(ent_array[loavgind:, 1] - avg)
y = mlab.normpdf(bins, mu, sigma)
l = plt.plot(bins, y, 'r--')
mu_magnitude = np.floor(np.log10(np.abs(mu)))
mu /= np.power(10, mu_magnitude)
sigma_magnitude = np.floor(np.log10(sigma))
sigma /= np.power(10, sigma_magnitude)
plt.figtext(0.965, 0.80,
            '$\mu = %.2f \cdot 10^{%i}$\n$\sigma = %.2f \cdot 10^{%i}$' % (mu, mu_magnitude, sigma, sigma_magnitude),
            ha='right', va='bottom', multialignment="left")
plt.xlabel(r'$\Delta S_{sub}$')
plt.ylabel(r'PD')
plt.tight_layout(padding)
plt.savefig(pltfolder + 'entropy_subsystem_fluctuations.eps', format='eps', dpi=1000)
plt.clf()
print('.', end='', flush=True)

print(" done!")


text:
import astropy.units as u
from astropy.units import UnitsError
from iminuit import Minuit
from iminuit.cost import UnbinnedNLL, BinnedNLL
import numpy as np
import pytest
from emmanoulopoulos.emmanoulopoulos_lc_simulation import power_spectral_density
from scipy.stats import norm, poisson, lognorm, gamma

def test_power_spectral_density():
    from emmanoulopoulos.models import power_spectral_density
    A = 0.0004325697931186756
    alpha_low = 1.7180692792824444
    alpha_high = 6.286876758053707
    f_bend = 0.014238670906664863
    c = 0.0020197057540771723

    f = 1e-3

    p = power_spectral_density(f, A, alpha_low, alpha_high, f_bend, c)

    np.testing.assert_array_almost_equal(p, 61.698685461074916, decimal=20)

def test_create_lc_uneven(uneven_times, sample_fluxes_from_bending_PL, tbin):
    from emmanoulopoulos.lightcurve import LC

    lc_original = LC(original_time=uneven_times, original_flux=sample_fluxes_from_bending_PL, errors=0.1*sample_fluxes_from_bending_PL, tbin=tbin)

    assert lc_original.interp_length == int((uneven_times.max() - uneven_times.min()) / tbin) + 1
    assert lc_original.interp_length == 205  # results from drawn random number with fixed known seed!

def test_create_lc_check_units_missing_unit(uneven_times, sample_fluxes_from_bending_PL):
    from emmanoulopoulos.lightcurve import LC
    tbin = 2
    with pytest.raises(TypeError) as err:
        lc_original = LC(original_time=uneven_times, original_flux=sample_fluxes_from_bending_PL, errors=0.1*sample_fluxes_from_bending_PL, tbin=tbin)
    assert err.type is TypeError

def test_create_lc_check_units_wrong_unit(uneven_times, sample_fluxes_from_bending_PL):
    from emmanoulopoulos.lightcurve import LC
    tbin = 2 * u.m
    with pytest.raises(UnitsError) as err:
        lc_original = LC(original_time=uneven_times, original_flux=sample_fluxes_from_bending_PL, errors=0.1*sample_fluxes_from_bending_PL, tbin=tbin)
    assert err.type is UnitsError

def test_create_lc_check_units_mismatching_units(uneven_times, sample_fluxes_from_bending_PL):
    from emmanoulopoulos.lightcurve import LC
    tbin = 2 * u.s
    with pytest.raises(ValueError) as err:
        lc_original = LC(original_time=uneven_times, original_flux=sample_fluxes_from_bending_PL, errors=0.1*sample_fluxes_from_bending_PL, tbin=tbin)
    assert err.type is ValueError

def test_create_lc_bad_length(uneven_times, sample_fluxes_from_bending_PL, tbin):
    from emmanoulopoulos.lightcurve import LC

    with pytest.raises(ValueError) as err:
        lc = LC(original_time=uneven_times[:-5], original_flux=sample_fluxes_from_bending_PL, errors=0.1*sample_fluxes_from_bending_PL, tbin=tbin)
    assert err.type is ValueError

def test_f_j(lc):
    f_j = lc.f_j()
    print(f_j)
    assert f_j[0] == 0
    assert len(f_j) == (lc.interp_length - lc.interp_length % 2) / 2 + 1

def test_lc_periodogram(lc):
    f_j, P_j = lc.periodogram()

    assert len(f_j) == len(P_j)
    assert len(f_j) == (lc.interp_length - lc.interp_length % 2) / 2 + 1

def test_fit_PSD(lc, psd_parameter):
    psd_parameter_fit = lc.fit_PSD()

    assert psd_parameter_fit["A"] is not None
    assert psd_parameter_fit["alpha_high"] is not None

def test_unbinned_fit_PDF(lc):
    pdf_unbinnned_fit = lc.fit_PDF(unbinned=True)

    assert pdf_unbinnned_fit.to_dict()["a"] > 0
    assert (pdf_unbinnned_fit.to_dict()["p"] > 0) and (pdf_unbinnned_fit.to_dict()["p"] < 1)

def test_binned_fit_PDF(lc):
    pdf_binnned_fit = lc.fit_PDF(unbinned=False)

    assert pdf_binnned_fit.to_dict()["a"] > 0
    assert (pdf_binnned_fit.to_dict()["p"] > 0) and (pdf_binnned_fit.to_dict()["p"] < 1)

max_stars_repo_name: lena-lin/emmanoulopoulos

max_issues_repo_name: lena-lin/emmanoulopoulos

max_forks_repo_name: lena-lin/emmanoulopoulos

reason: import numpy,from scipy,import astropy,from astropy


text:
import numpy as np
import cv2
import matplotlib.pyplot as plt
from MovingObjectDetector.BackgroundModel import BackgroundModel
from MovingObjectDetector.DetectionRefinement import DetectionRefinement
import TrainNetwork.BaseFunctions as basefunctions

input_image_idx = 277

image_idx_offset = 8026480
num_of_template = 4
imagefolder = "E:/WASABI-AngelFire-02/"

image_idx_offset = 0
num_of_template = 4

model_binary, aveImg_binary, model_regression, aveImg_regression = basefunctions.ReadModels()

images = []
for i in range(num_of_template):
    ReadImage = cv2.imread(imagefolder + "frame%06d.png" % (input_image_idx+image_idx_offset+i-num_of_template), cv2.IMREAD_GRAYSCALE)
    images.append(ReadImage)
bgt = BackgroundModel(num_of_template=num_of_template, templates=images)
input_image = cv2.imread(imagefolder + "frame%06d.png" % (input_image_idx+image_idx_offset), cv2.IMREAD_GRAYSCALE)
Hs = bgt.doCalculateHomography(input_image)
bgt.doMotionCompensation(Hs, input_image.shape)
BackgroundSubtractionCentres, BackgroundSubtractionProperties = bgt.doBackgroundSubtraction(input_image, thres=10)

dr = DetectionRefinement(input_image, bgt.getCompensatedImages(), BackgroundSubtractionCentres, BackgroundSubtractionProperties, model_binary, aveImg_binary, model_regression, aveImg_regression)
refinedDetections, refinedProperties = dr.doMovingVehicleRefinement()
regressedDetections = dr.doMovingVehiclePositionRegression()
regressedDetections = np.asarray(regressedDetections)

plt.figure(11)
plt.imshow(np.repeat(np.expand_dims(input_image, -1), 3, axis=2))
plt.plot(BackgroundSubtractionCentres[:,0], BackgroundSubtractionCentres[:,1], 'g.')
plt.plot(refinedDetections[:,0], refinedDetections[:,1], 'y.')
plt.plot(np.int32(regressedDetections[:,0]), np.int32(regressedDetections[:,1]), 'r.')
plt.show()


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

import numpy as np
import tensorflow as tf

import logging

logger = logging.getLogger(__name__)

def clip_sentence(sentence, sizes):
    """
    Clip the input sentence placeholders to the length of the longest one in the batch.
    This saves processing time.

    :param sentence: tensor with shape (batch, time_steps)
    :param sizes: tensor with shape (batch)
    :return: tensor with shape (batch, time_steps)
    """
    return tf.slice(sentence, [0, 0], tf.stack([-1, tf.reduce_max(sizes)]))

def attention_softmax3d(values):
    """
    Performs a softmax over the attention values.
    :param values: tensor with shape (batch_size, time_steps, time_steps)
    :return: tensor with shape (batch_size, time_steps, time_steps)
    """
    original_shape = tf.shape(values)
    # tensor with shape (batch_size * time_steps, time_steps)
    reshaped_values = tf.reshape(tensor=values, shape=[-1, original_shape[2]])
    # tensor with shape (batch_size * time_steps, time_steps)
    softmax_reshaped_values = tf.nn.softmax(reshaped_values)
    # tensor with shape (batch_size, time_steps, time_steps)
    return tf.reshape(softmax_reshaped_values, original_shape)

def mask_3d(sequences, sequence_lengths, mask_value, dimension=2):
    """
    Given a batch of matrices, each with shape m x n, mask the values in each
    row after the positions indicated in sentence_sizes.
    This function is supposed to mask the last columns in the raw attention
    matrix (e_{i, j}) in cases where the sentence2 is smaller than the
    maximum.
    :param sequences: tensor with shape (batch_size, m, n)
    :param sequence_lengths: tensor with shape (batch_size) containing the sentence sizes that
        should be limited
    :param mask_value: scalar value to assign to items after sentence size
    :param dimension: over which dimension to mask values
    :return: a tensor with the same shape as `values`
    """
    if dimension == 1:
        sequences = tf.transpose(sequences, [0, 2, 1])
    time_steps1, time_steps2 = tf.shape(sequences)[1], tf.shape(sequences)[2]
    ones = tf.ones_like(sequences, dtype=tf.int32)
    pad_values = mask_value * tf.cast(ones, tf.float32)
    mask = tf.sequence_mask(sequence_lengths, time_steps2)
    # mask is (batch_size, sentence2_size). we have to tile it for 3d
    mask3d = tf.tile(tf.expand_dims(mask, 1), (1, time_steps1, 1))
    masked = tf.where(mask3d, sequences, pad_values)
    return tf.transpose(masked, [0, 2, 1]) if dimension == 1 else masked

def distance_biases(time_steps, window_size=10, reuse=False):
    """
    Return a 2-d tensor with the values of the distance biases to be applied
    on the intra-attention matrix of size sentence_size

    :param time_steps: tensor scalar
    :param window_size: window size
    :param reuse: reuse variables
    :return: 2-d tensor (time_steps, time_steps)
    """
    with tf.variable_scope('distance-bias', reuse=reuse):
        # this is d_{i-j}
        distance_bias = tf.get_variable('dist_bias', [window_size], initializer=tf.zeros_initializer())
        r = tf.range(0, time_steps)
        r_matrix = tf.tile(tf.reshape(r, [1, -1]), tf.stack([time_steps, 1]))
        raw_idxs = r_matrix - tf.reshape(r, [-1, 1])
        clipped_idxs = tf.clip_by_value(raw_idxs, 0, window_size - 1)
        values = tf.nn.embedding_lookup(distance_bias, clipped_idxs)
    return values

def intra_attention(sequence, reuse=False):
    """
    Compute the intra attention of a sentence. It returns a concatenation
    of the original sentence with its attended output.

    :param sequence: tensor in shape (batch, time_steps, num_units)
    :param reuse: reuse variables
    :return: a tensor in shape (batch, time_steps, 2*num_units)
    """
    time_steps = tf.shape(sequence)[1]
    with tf.variable_scope('intra-attention') as _:
        # this is F_intra in the paper
        # f_intra1 is (batch, time_steps, num_units) and
        # f_intra1_t is (batch, num_units, time_steps)
        f_intra_t = tf.transpose(sequence, [0, 2, 1])

        # these are f_ij
        # raw_attentions is (batch, time_steps, time_steps)
        raw_attentions = tf.matmul(sequence, f_intra_t)

        # bias has shape (time_steps, time_steps)
        bias = distance_biases(time_steps, reuse=reuse)

        # bias is broadcast along batches
        raw_attentions += bias
        attentions = attention_softmax3d(raw_attentions)
        attended = tf.matmul(attentions, sequence)
    return tf.concat(axis=2, values=[sequence, attended])

def count_trainable_parameters(var_list=None):
    """
    Count the number of trainable tensorflow parameters loaded in
    the current graph.
    """
    total_params = 0
    if var_list is None:
        var_list = tf.trainable_variables()
    for variable in var_list:
        variable_params = np.prod([1] + [dim.value for dim in variable.get_shape()])
        logging.info('{0} ({1}): {2} params'.format(variable.name, str(variable.get_shape()), variable_params))
        total_params += variable_params
    return total_params

def get_variables_in_scope(scope_name):
    return tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope=scope_name)

max_stars_repo_name: issca/inferbeddings

max_issues_repo_name: issca/inferbeddings

max_forks_repo_name: issca/inferbeddings


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# LICENSE
#
# Copyright (C) 2010-2018 GEM Foundation, G. Weatherill, M. Pagani,
# D. Monelli.
#
# The Hazard Modeller's Toolkit is free software: you can redistribute
# it and/or modify it under the terms of the GNU Affero General Public
# License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>
#
# DISCLAIMER
#
# The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein
# is released as a prototype implementation on behalf of
# scientists and engineers working within the GEM Foundation (Global
# Earthquake Model).
#
# It is distributed for the purpose of open collaboration and in the
# hope that it will be useful to the scientific, engineering, disaster
# risk and software design communities.
#
# The software is NOT distributed as part of GEM's OpenQuake suite
# (https://www.globalquakemodel.org/tools-products) and must be considered as a
# separate entity. The software provided herein is designed and implemented
# by scientific staff. It is not developed to the design standards, nor
# subject to same level of critical review by professional software
# developers, as GEM's OpenQuake software suite.
#
# Feedback and contribution to the software is welcome, and can be
# directed to the hazard scientific staff of the GEM Model Facility
# (hazard@globalquakemodel.org).
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# The GEM Foundation, and the authors of the software, assume no
# liability for use of the software.

# -*- coding: utf-8 -*-

'''
Tests the construction and methods within the :class:
openquake.hmtk.sources.area_source.mtkAreaSource
'''

import unittest
import warnings
import numpy as np
from openquake.hazardlib.geo import point, polygon
from openquake.hazardlib.source.area import AreaSource
from openquake.hazardlib.pmf import PMF
from openquake.hazardlib.tom import PoissonTOM
from openquake.hazardlib.mfd.truncated_gr import TruncatedGRMFD
from openquake.hazardlib.scalerel.wc1994 import WC1994
from openquake.hmtk.sources.area_source import mtkAreaSource
from openquake.hmtk.seismicity.catalogue import Catalogue
from openquake.hmtk.seismicity.selector import CatalogueSelector

TOM = PoissonTOM(50.0)
SOURCE_ATTRIBUTES = ['mfd', 'name', 'geometry', 'nodal_plane_dist', 'typology',
                     'upper_depth', 'catalogue', 'rupt_aspect_ratio',
                     'lower_depth', 'id', 'hypo_depth_dist', 'mag_scale_rel',
                     'trt']

class TestAreaSource(unittest.TestCase):
    '''
    Tester class for openquake.hmtk.sources.area_source.mtkAreaSource
    '''

    def setUp(self):
        warnings.simplefilter("ignore")
        self.catalogue = Catalogue()
        self.area_source = mtkAreaSource('101', 'A Source')

    def test_area_source_instantiation(self):
        # Tests the core (minimal) instantiation of the class
        # Check source has all required attributes
        self.assertListEqual(sorted(self.area_source.__dict__),
                             sorted(SOURCE_ATTRIBUTES))
        self.assertEqual(self.area_source.id, '101')
        self.assertEqual(self.area_source.name, 'A Source')
        self.assertEqual(self.area_source.typology, 'Area')

    def test_depth_checker(self):
        # Tests the checker to ensure correct depth values
        # Bad Case - Negative upper depths
        with self.assertRaises(ValueError) as ver:
            self.area_source._check_seismogenic_depths(-1.0, 20.)
        self.assertEqual(str(ver.exception),
                         'Upper seismogenic depth must be greater than or '
                         'equal to 0.0!')

        # Bad Case - Lower depth smaller than upper depth
        with self.assertRaises(ValueError) as ver:
            self.area_source._check_seismogenic_depths(30., 20.)
        self.assertEqual(str(ver.exception),
                         'Lower seismogenic depth must take a greater value '
                         'than upper seismogenic depth')
        # Good Case
        self.area_source._check_seismogenic_depths(0.0, 20.)
        self.assertAlmostEqual(0.0, self.area_source.upper_depth)
        self.assertAlmostEqual(20.0, self.area_source.lower_depth)

    def test_geometry_inputs(self):
        # Tests the geometry definition
        simple_polygon = polygon.Polygon([point.Point(2.0, 3.0),
                                          point.Point(3.0, 3.0),
                                          point.Point(3.0, 2.0),
                                          point.Point(2.0, 2.0)])

        simple_polygon_array = np.array([[2.0, 3.0],
                                         [3.0, 3.0],
                                         [3.0, 2.0],
                                         [2.0, 2.0]])
        # Using nhlib.geo.polygon.Polygon class as input
        self.area_source.create_geometry(simple_polygon, 0.0, 30.0)
        # Check that geometry is an instance of nhlib.geo.polygon.Polygon
        self.assertTrue(isinstance(self.area_source.geometry,
                                   polygon.Polygon))

        np.testing.assert_array_almost_equal(self.area_source.geometry.lons,
                                             np.array([2., 3., 3., 2.]))
        np.testing.assert_array_almost_equal(self.area_source.geometry.lats,
                                             np.array([3., 3., 2., 2.]))
        self.assertAlmostEqual(0.0, self.area_source.upper_depth)
        self.assertAlmostEqual(30.0, self.area_source.lower_depth)

        self.area_source = mtkAreaSource('101', 'A Source')
        # Using numpy array as input
        self.area_source.create_geometry(simple_polygon_array, 0.0, 30.0)
        self.assertTrue(isinstance(self.area_source.geometry,
                                   polygon.Polygon))

        # Check that geometry is an instance of nhlib.geo.polygon.Polygon
        np.testing.assert_array_almost_equal(self.area_source.geometry.lons,
                                             np.array([2., 3., 3., 2.]))
        np.testing.assert_array_almost_equal(self.area_source.geometry.lats,
                                             np.array([3., 3., 2., 2.]))

        self.assertAlmostEqual(0.0, self.area_source.upper_depth)
        self.assertAlmostEqual(30.0, self.area_source.lower_depth)

        # For any other input type - check ValueError is raised
        self.area_source = mtkAreaSource('101', 'A Source')
        with self.assertRaises(ValueError) as ver:
            self.area_source.create_geometry('a bad input', 0.0, 30.0)
        self.assertEqual(str(ver.exception),
                         'Unrecognised or unsupported geometry definition')

        # For numpy array with only two rows
        self.area_source = mtkAreaSource('101', 'A Source')
        simple_polygon_array = np.array([[2.0, 3.0],
                                         [3.0, 3.0]])
        with self.assertRaises(ValueError) as ver:
            self.area_source.create_geometry(simple_polygon_array, 0.0, 30.0)
        self.assertEqual(str(ver.exception),
                         'Incorrectly formatted polygon geometry - '
                         'needs three or more vertices')

    def test_select_events_in_source(self):
        # Basic test of method to select events from catalogue in polygon
        self.area_source = mtkAreaSource('101', 'A Source')
        simple_polygon = polygon.Polygon([point.Point(2.0, 3.0),
                                          point.Point(3.0, 3.0),
                                          point.Point(3.0, 2.0),
                                          point.Point(2.0, 2.0)])

        self.catalogue.data['eventID'] = np.arange(0, 7, 1)
        self.catalogue.data['longitude'] = np.arange(1.0, 4.5, 0.5)
        self.catalogue.data['latitude'] = np.arange(1.0, 4.5, 0.5)
        self.catalogue.data['depth'] = np.ones(7, dtype=float)
        # Simple Case - No buffer
        selector0 = CatalogueSelector(self.catalogue)
        self.area_source.create_geometry(simple_polygon, 0., 30.)
        self.area_source.select_catalogue(selector0, 0.)
        np.testing.assert_array_almost_equal(
            np.array([2., 2.5, 3.]),
            self.area_source.catalogue.data['longitude'])

        np.testing.assert_array_almost_equal(
            np.array([2., 2.5, 3.]),
            self.area_source.catalogue.data['latitude'])

        np.testing.assert_array_almost_equal(
            np.array([1., 1., 1.]),
            self.area_source.catalogue.data['depth'])

        # Simple case - dilated by 200 km (selects all)
        self.area_source.select_catalogue(selector0, 200.)
        np.testing.assert_array_almost_equal(
            np.array([1., 1.5, 2., 2.5, 3., 3.5, 4.0]),
            self.area_source.catalogue.data['longitude'])

        np.testing.assert_array_almost_equal(
            np.array([1., 1.5, 2., 2.5, 3., 3.5, 4.0]),
            self.area_source.catalogue.data['latitude'])

        np.testing.assert_array_almost_equal(
            np.ones(7, dtype=float),
            self.area_source.catalogue.data['depth'])

        # Bad case - no events in catalogue
        self.catalogue = Catalogue()
        selector0 = CatalogueSelector(self.catalogue)

        with self.assertRaises(ValueError) as ver:
            self.area_source.select_catalogue(selector0, 0.0)
            self.assertEqual(str(ver.exception),
                             'No events found in catalogue!')

    def test_create_oqhazardlib_source(self):
        # Define a complete source
        area_geom = polygon.Polygon([point.Point(10., 10.),
                                     point.Point(12., 10.),
                                     point.Point(12., 8.),
                                     point.Point(10., 8.)])
        mfd1 = TruncatedGRMFD(5.0, 8.0, 0.1, 3.0, 1.0)
        self.area_source = mtkAreaSource(
            '001',
            'A Point Source',
            trt='Active Shallow Crust',
            geometry=area_geom,
            upper_depth=0.,
            lower_depth=20.,
            mag_scale_rel=None,
            rupt_aspect_ratio=1.0,
            mfd=mfd1,
            nodal_plane_dist=None,
            hypo_depth_dist=None)
        test_source = self.area_source.create_oqhazardlib_source(TOM, 1.0,
                                                                 10.0, True)
        self.assertIsInstance(test_source, AreaSource)
        self.assertIsInstance(test_source.mfd, TruncatedGRMFD)
        self.assertAlmostEqual(test_source.mfd.b_val, 1.0)
        self.assertIsInstance(test_source.nodal_plane_distribution, PMF)
        self.assertIsInstance(test_source.hypocenter_distribution, PMF)
        self.assertIsInstance(test_source.magnitude_scaling_relationship,
                              WC1994)

    def tearDown(self):
        warnings.resetwarnings()

max_stars_repo_name: gfzriesgos/shakyground-lfs

max_issues_repo_name: gfzriesgos/shakyground-lfs

max_forks_repo_name: gfzriesgos/shakyground-lfs