# Path: pactools/utils/maths.py
# def square(c):
#     """square a complex array"""
#     return np.real(c * np.conjugate(c))
# 
# def is_power2(num):
#     """Test if num is a power of 2. (int -> bool)"""
#     num = int(num)
#     return num != 0 and ((num & (num - 1)) == 0)
# 
# def prime_factors(num):
#     """
#     Decomposition in prime factor.
#     Used to find signal length that speed up Hilbert transform.
# 
#     Parameters
#     ----------
#     num : int
# 
#     Returns
#     -------
#     decomposition : list of int
#         List of prime factors sorted by ascending order
#     """
#     assert num > 0
#     assert isinstance(num, int)
# 
#     decomposition = []
#     for i in chain([2], range(3, int(np.sqrt(num)) + 1, 2)):
#         while num % i == 0:
#             num = num // i
#             decomposition.append(i)
#         if num == 1:
#             break
#     if num != 1:
#         decomposition.append(num)
#     return decomposition
# 
# def compute_n_fft(signals):
#     """
#     Compute the number of element in the FFT to have a good prime
#     decomposition, for hilbert filter.
#     Parameters
#     ----------
#     signals : array
# 
#     Returns
#     -------
#     n_fft: integer >= signals.shape[-1]
#     """
#     n_fft = signals.shape[-1]
#     while prime_factors(n_fft)[-1] > 20:
#         n_fft += 1
#     return n_fft
# 
# def next_power2(num):
#     """Compute the smallest power of 2 >= to num.(float -> int)"""
#     return 2 ** int(np.ceil(np.log2(num)))
# 
# Path: pactools/utils/viz.py
# def compute_vmin_vmax(spec, vmin=None, vmax=None, tick=0.01, percentile=1):
#     """compute automatic scale for plotting `spec`"""
#     if percentile > 100:
#         percentile = 100
#     if percentile < 0:
#         percentile = 0
# 
#     if vmin is None:
#         vmin = np.floor(np.percentile(spec, percentile) / tick) * tick
#     if vmax is None:
#         vmax = np.ceil(np.percentile(spec, 100 - percentile) / tick) * tick
# 
#     return vmin, vmax


import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.linalg import hankel
from scipy.signal import hilbert
from .maths import square, is_power2, prime_factors, compute_n_fft, next_power2
from .viz import compute_vmin_vmax
        bicoherence[np.triu_indices(n_freq, 1)] = 0
        bicoherence = bicoherence[:, :n_freq // 2 + 1]

        vmin, vmax = compute_vmin_vmax(bicoherence, tick=1e-15, percentile=1)
        ax.imshow(bicoherence, cmap=plt.cm.viridis, aspect='auto', vmin=vmin,
                  vmax=vmax, origin='lower', extent=(0, fmax // 2, 0, fmax),
                  interpolation='none')
        ax.set_title('Bicoherence (%s)' % self.method)
        ax.set_xlabel('Frequency (Hz)')
        ax.set_ylabel('Frequency (Hz)')
        # add_colorbar(fig, cax, vmin, vmax, unit='', ax=ax)
        return ax

    def main_frequency(self):
        pass


def phase_amplitude(signals, phase=True, amplitude=True):
    """Extract instantaneous phase and amplitude with Hilbert transform"""
    # one dimension array
    if signals.ndim == 1:
        signals = signals[None, :]
        one_dim = True
    elif signals.ndim == 2:
        one_dim = False
    else:
        raise ValueError('Impossible to compute phase_amplitude with ndim ='
                         ' %s.' % (signals.ndim, ))

    n_epochs, n_points = signals.shape

# Path: pactools/utils/maths.py
# def square(c):
#     """square a complex array"""
#     return np.real(c * np.conjugate(c))
# 
# def is_power2(num):
#     """Test if num is a power of 2. (int -> bool)"""
#     num = int(num)
#     return num != 0 and ((num & (num - 1)) == 0)
# 
# def prime_factors(num):
#     """
#     Decomposition in prime factor.
#     Used to find signal length that speed up Hilbert transform.
# 
#     Parameters
#     ----------
#     num : int
# 
#     Returns
#     -------
#     decomposition : list of int
#         List of prime factors sorted by ascending order
#     """
#     assert num > 0
#     assert isinstance(num, int)
# 
#     decomposition = []
#     for i in chain([2], range(3, int(np.sqrt(num)) + 1, 2)):
#         while num % i == 0:
#             num = num // i
#             decomposition.append(i)
#         if num == 1:
#             break
#     if num != 1:
#         decomposition.append(num)
#     return decomposition
# 
# def compute_n_fft(signals):
#     """
#     Compute the number of element in the FFT to have a good prime
#     decomposition, for hilbert filter.
#     Parameters
#     ----------
#     signals : array
# 
#     Returns
#     -------
#     n_fft: integer >= signals.shape[-1]
#     """
#     n_fft = signals.shape[-1]
#     while prime_factors(n_fft)[-1] > 20:
#         n_fft += 1
#     return n_fft
# 
# def next_power2(num):
#     """Compute the smallest power of 2 >= to num.(float -> int)"""
#     return 2 ** int(np.ceil(np.log2(num)))
# 
# Path: pactools/utils/viz.py
# def compute_vmin_vmax(spec, vmin=None, vmax=None, tick=0.01, percentile=1):
#     """compute automatic scale for plotting `spec`"""
#     if percentile > 100:
#         percentile = 100
#     if percentile < 0:
#         percentile = 0
# 
#     if vmin is None:
#         vmin = np.floor(np.percentile(spec, percentile) / tick) * tick
#     if vmax is None:
#         vmax = np.ceil(np.percentile(spec, 100 - percentile) / tick) * tick
# 
#     return vmin, vmax


# Path: pactools/utils/spectrum.py
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.linalg import hankel
from scipy.signal import hilbert
from .maths import square, is_power2, prime_factors, compute_n_fft, next_power2
from .viz import compute_vmin_vmax

        bicoherence[np.triu_indices(n_freq, 1)] = 0
        bicoherence = bicoherence[:, :n_freq // 2 + 1]

        vmin, vmax = compute_vmin_vmax(bicoherence, tick=1e-15, percentile=1)
        ax.imshow(bicoherence, cmap=plt.cm.viridis, aspect='auto', vmin=vmin,
                  vmax=vmax, origin='lower', extent=(0, fmax // 2, 0, fmax),
                  interpolation='none')
        ax.set_title('Bicoherence (%s)' % self.method)
        ax.set_xlabel('Frequency (Hz)')
        ax.set_ylabel('Frequency (Hz)')
        # add_colorbar(fig, cax, vmin, vmax, unit='', ax=ax)
        return ax

    def main_frequency(self):
        pass


def phase_amplitude(signals, phase=True, amplitude=True):
    """Extract instantaneous phase and amplitude with Hilbert transform"""
    # one dimension array
    if signals.ndim == 1:
        signals = signals[None, :]
        one_dim = True
    elif signals.ndim == 2:
        one_dim = False
    else:
        raise ValueError('Impossible to compute phase_amplitude with ndim ='
                         ' %s.' % (signals.ndim, ))

    n_epochs, n_points = signals.shape

    n_fft = compute_n_fft(signals)