# Path: dppy/beta_ensemble_polynomial_potential_core.py
# def shift_pol(P, x0):
# def log_pdf_convex_quartic(x, P):
# def find_mode_convex_quartic(P):
# def find_b_a_convex_quartic(P, mode):
# def log_pdf_convex_gen_gamma(x, shape, P):
# def find_mode_convex_gen_gamma(shape, P):
# def find_a_b_convex_gen_gamma(shape, P, mode):
#     def f(x):
#     def d1_f(x):
#     def d2_f(x):
# def trunc_gaussian_sampler(b, mu, var, random_state=None):
# def gen_gamma_alpha_lt_1_sampler(shape, P, random_state=None):
# def sampler_exact_convex_quartic(P, shape=None, random_state=None):
#         def log_f(x):
#         def log_f(x):
#     def log_h(x):
# def sampler_mala(x, V, sigma=0.01, nb_steps=20, random_state=None):
# def polynomial_in_negative_log_conditional_a_coef(i, a, b, V):
# def polynomial_in_negative_log_conditional_b_coef(i, a, b, V):
# def equilibrium_x2_x4(g_2, g_4):
#         def dmu(x):
#         def dmu(x):
# def cdf_equilibrium_x2_x4(g_2, g_4):
#     def cdf(x):
# def equilibrium_x2m(m, g_2m):
#     def mu_eq(x):
# def cdf_equilibrium_x2m(m, g_2m):
#         def cdf(x):
#         def cdf(x):
# def equilibrium_ClItKr10():
#     def mu_eq(x):
#     def __init__(self):
#     def change_of_variable(self, s):  # , old_implem=False):
#         def phi(x):
#         def d_phi(x):
#     def kernel(self, s):  # , old_implem=False):
#         def K(X, Y=None, eps_xy=1e-5):
#     def cdf(self, x, quad_order=50):  # , old_implem=False):
#     def airy_kernel(x, y, eps_xy=1e-5):
#     def compute_quadrature(order, type='Legendre'):
#     def fredholm_determinant(kernel, x_quad, w_quad):  # , old_implem):
#         X = rng.normal(loc=b, scale=np.sqrt(var))
#         U = rng.rand()
#         X = 1e-7
#         X = rng.gamma(shape=shape, scale=1 / P[1], size=1)
#         U = rng.rand()
#         U, V, W = rng.rand(3)  # Unif_[0,1]
#             X = _b + (_a - _b) * U
#             X = _a + _a * U
#             X = _b + _b * U
#             X = _2a - mu * np.log(U)
#             X = _2b + mu * np.log(1.0 - U * (1.0 - np.exp(-trunc / mu)))
#     P = np.poly1d(0.0)
#     P = np.poly1d(0.0)
#         X = np.arcsin(x / a)
#             X = np.arcsin(x / a)
#             X = np.arcsin(x / a)
#             X = np.ravel(X)
#         K = kernel(x_quad, x_quad)
# class TracyWidom(object):


import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
from numpy.polynomial.chebyshev import poly2cheb  # cf rescale_largest_eig_val
from dppy.beta_ensemble_polynomial_potential_core import (polynomial_in_negative_log_conditional_a_coef as P_a_cond,
    polynomial_in_negative_log_conditional_b_coef as P_b_cond,
    sampler_exact_convex_quartic,
    sampler_mala,
    equilibrium_x2m,
    equilibrium_x2_x4,
    equilibrium_ClItKr10,
    check_random_state)
        self.nb_gibbs_passes = nb_gibbs_passes

        a, b = np.zeros((2, N + 3))

        if return_chain_of_eig_vals:
            eig_vals = np.zeros((N, nb_gibbs_passes))
        elif return_chain_of_lambda_max:
            lambda_max = np.zeros(nb_gibbs_passes)

        for p in range(nb_gibbs_passes):
            if (p + 1) % 50 == 0:
                print(p + 1)

            for i in range(1, N + 1):

                # a_i | ... propto exp - P_a_i
                P_a_i = 0.5 * self.beta * N * P_a_cond(i, a, b, self.V)
                if sample_exact_cond:
                    a[i], _ = sampler_exact_convex_quartic(
                                P=P_a_i,
                                random_state=rng)
                else:
                    a[i] = sampler_mala(a[i],
                                        V=P_a_i,
                                        sigma=0.01,
                                        nb_steps=nb_mala_steps,
                                        random_state=rng)

                # b_i | ... propto x^(shape-1) * exp - P_b_i
                if i < N:

# Path: dppy/beta_ensemble_polynomial_potential_core.py
# def shift_pol(P, x0):
# def log_pdf_convex_quartic(x, P):
# def find_mode_convex_quartic(P):
# def find_b_a_convex_quartic(P, mode):
# def log_pdf_convex_gen_gamma(x, shape, P):
# def find_mode_convex_gen_gamma(shape, P):
# def find_a_b_convex_gen_gamma(shape, P, mode):
#     def f(x):
#     def d1_f(x):
#     def d2_f(x):
# def trunc_gaussian_sampler(b, mu, var, random_state=None):
# def gen_gamma_alpha_lt_1_sampler(shape, P, random_state=None):
# def sampler_exact_convex_quartic(P, shape=None, random_state=None):
#         def log_f(x):
#         def log_f(x):
#     def log_h(x):
# def sampler_mala(x, V, sigma=0.01, nb_steps=20, random_state=None):
# def polynomial_in_negative_log_conditional_a_coef(i, a, b, V):
# def polynomial_in_negative_log_conditional_b_coef(i, a, b, V):
# def equilibrium_x2_x4(g_2, g_4):
#         def dmu(x):
#         def dmu(x):
# def cdf_equilibrium_x2_x4(g_2, g_4):
#     def cdf(x):
# def equilibrium_x2m(m, g_2m):
#     def mu_eq(x):
# def cdf_equilibrium_x2m(m, g_2m):
#         def cdf(x):
#         def cdf(x):
# def equilibrium_ClItKr10():
#     def mu_eq(x):
#     def __init__(self):
#     def change_of_variable(self, s):  # , old_implem=False):
#         def phi(x):
#         def d_phi(x):
#     def kernel(self, s):  # , old_implem=False):
#         def K(X, Y=None, eps_xy=1e-5):
#     def cdf(self, x, quad_order=50):  # , old_implem=False):
#     def airy_kernel(x, y, eps_xy=1e-5):
#     def compute_quadrature(order, type='Legendre'):
#     def fredholm_determinant(kernel, x_quad, w_quad):  # , old_implem):
#         X = rng.normal(loc=b, scale=np.sqrt(var))
#         U = rng.rand()
#         X = 1e-7
#         X = rng.gamma(shape=shape, scale=1 / P[1], size=1)
#         U = rng.rand()
#         U, V, W = rng.rand(3)  # Unif_[0,1]
#             X = _b + (_a - _b) * U
#             X = _a + _a * U
#             X = _b + _b * U
#             X = _2a - mu * np.log(U)
#             X = _2b + mu * np.log(1.0 - U * (1.0 - np.exp(-trunc / mu)))
#     P = np.poly1d(0.0)
#     P = np.poly1d(0.0)
#         X = np.arcsin(x / a)
#             X = np.arcsin(x / a)
#             X = np.arcsin(x / a)
#             X = np.ravel(X)
#         K = kernel(x_quad, x_quad)
# class TracyWidom(object):


# Path: dppy/beta_ensemble_polynomial_potential.py
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
from numpy.polynomial.chebyshev import poly2cheb  # cf rescale_largest_eig_val
from dppy.beta_ensemble_polynomial_potential_core import (polynomial_in_negative_log_conditional_a_coef as P_a_cond,
    polynomial_in_negative_log_conditional_b_coef as P_b_cond,
    sampler_exact_convex_quartic,
    sampler_mala,
    equilibrium_x2m,
    equilibrium_x2_x4,
    equilibrium_ClItKr10,
    check_random_state)

        self.nb_gibbs_passes = nb_gibbs_passes

        a, b = np.zeros((2, N + 3))

        if return_chain_of_eig_vals:
            eig_vals = np.zeros((N, nb_gibbs_passes))
        elif return_chain_of_lambda_max:
            lambda_max = np.zeros(nb_gibbs_passes)

        for p in range(nb_gibbs_passes):
            if (p + 1) % 50 == 0:
                print(p + 1)

            for i in range(1, N + 1):

                # a_i | ... propto exp - P_a_i
                P_a_i = 0.5 * self.beta * N * P_a_cond(i, a, b, self.V)
                if sample_exact_cond:
                    a[i], _ = sampler_exact_convex_quartic(
                                P=P_a_i,
                                random_state=rng)
                else:
                    a[i] = sampler_mala(a[i],
                                        V=P_a_i,
                                        sigma=0.01,
                                        nb_steps=nb_mala_steps,
                                        random_state=rng)

                # b_i | ... propto x^(shape-1) * exp - P_b_i
                if i < N:

                    P_b_i = 0.5 * self.beta * N * P_b_cond(i, a, b, self.V)