def spectrum(x,
             n0=0,
             T_s=1,
             oversample=1,
             only_positive=True):
    """
    Return the spectrum for the signal *x* calculated via FFT and the
    associated frequencies as a tuple. The *n0* parameter gives the
    index in *x* for time index 0 (*n0* = 0 means that `x[0]` is at
    time 0). The number of spectral samples returned is the next power
    of 2 greater than the length of *x* multiplied by *oversample*. If
    *only_positive*,return the spectrum only for positive frequencies
    (and raise an exception if *x* is not real).
    """
    assert oversample >= 1 and isinstance(oversample, int)
    N = nextpow2(len(x)) * 2**(oversample - 1)
    X = NP.fft.fft(x, n=N) * T_s
    f = NP.fft.fftfreq(N, d=T_s)
    if n0 != 0:
        X *= NP.exp(-1j * 2 * math.pi * NP.arange(N) * n0 / N)
    X = NP.fft.fftshift(X)
    f = NP.fft.fftshift(f)
    if only_positive:
        if any(NP.iscomplex(x)):
            raise ValueError(''x is complex and only returning information for positive frequencies --- this is likely not what you want to do'')
        X = X[f >= 0]
        f = f[f >= 0]
    return X, f

def plot_data(data, scroll_axis=2):
    """ Plot an image associated data.
    Currently support on 1D,2D or 3D data.
 
    Parameters
    ----------
    data: array
        the data to be displayed.
    scroll_axis: int (optional,default 2)
        the scroll axis for 3d data.
    """
    # Check input parameters
    if data.ndim not in range(1, 4):
        raise ValueError("Unsupported data dimension.")
 
    # Deal with complex data
    if numpy.iscomplex(data).any():
        data = numpy.abs(data)
 
    # Create application
    app = pyqtgraph.mkQApp()
 
    # Create the widget
    if data.ndim == 3:
        indices = [i for i in range(3) if i != scroll_axis]
        indices = [scroll_axis] + indices
        widget = pyqtgraph.image(numpy.transpose(data, indices))
    elif data.ndim == 2:
        widget = pyqtgraph.image(data)
    else:
        widget = pyqtgraph.plot(data)
 
    # Run application
    app.exec_()

def average_structure(X):
    """
    Calculate an average structure from an ensemble of structures
    (i.e. X is a rank-3 tensor: X[i] is a (N,3) configuration matrix).
 
    @param X: m x n x 3 input vector
    @type X: numpy array
 
    @return: average structure
    @rtype: (n,3) numpy.array
    """
    from numpy.linalg import eigh
 
    B = csb.numeric.gower_matrix(X)
    v, U = eigh(B)
    if numpy.iscomplex(v).any():
        v = v.real
    if numpy.iscomplex(U).any():
        U = U.real
 
    indices = numpy.argsort(v)[-3:]
    v = numpy.take(v, indices, 0)
    U = numpy.take(U, 1)
 
    x = U * numpy.sqrt(v)
    i = 0
    while is_mirror_image(x, X[0]) and i < 2:
        x[:, i] *= -1
        i += 1
    return x

def rel_entropy_true(p, q):
    """KL divergence (relative entropy) D(p||q) in bits
 
    Returns a scalar entropy when the input distributions p and q are
    vectors of probability masses,or returns in a row vector the
    columnwise relative entropies of the input probability matrices p and
    q"""
 
    if type(p) == list or type(q) == tuple:
        p = np.array(p)
    if type(q) == list or type(q) == tuple:
        q = np.array(q)
 
    if not p.shape == q.shape:
        raise Exception(''p and q must be equal sizes'',
                        ''p: '' + str(p.shape),
                        ''q: '' + str(q.shape))
 
    if (np.iscomplex(p).any() or not
        np.isfinite(p).any() or
        (p < 0).any() or
        (p > 1).any()):
        raise Exception(''The probability elements of p must be real numbers''
                        ''between 0 and 1.'')
 
    if (np.iscomplex(q).any() or not
        np.isfinite(q).any() or
        (q < 0).any() or
        (q > 1).any()):
        raise Exception(''The probability elements of q must be real numbers''
                        ''between 0 and 1.'')
 
    eps = math.sqrt(np.spacing(1))
    if (np.abs(np.sum(p, axis=0) - 1) > eps).any():
        raise Exception(''Sum of the probability elements of p must equal 1.'')
    if (np.abs(np.sum(q, axis=0) - 1) > eps).any():
        raise Exception(''Sum of the probability elements of q must equal 1.'')
 
    return sum(np.log2((p**p) / (q**p)))

def pnorm(self, p):
        r"""
        Calculates the p-norm of a vector.
 
        Parameters
        ----------
        p : int
            Used in computing the p-norm. This should only be set
            when calculating the pnorm of a vector is desired.
 
        Returns
        -------
        pn : float
            The p-norm of the vector.
 
        Notes
        -----
        The p-norm,which is considered a class of vector norms is defined as:
 
        .. math::
 
            ||x||_p = \\sqrt[p]{|x_1|^p + |x_2|^p + \\cdots + |x_n|^p} \\qquad p \\geq 1
 
        References
        ----------
        Golub,G.,& Van Loan,C. (2013). Matrix computations (3rd ed.). Baltimore (MD): Johns Hopkins U.P.
 
        """
        if p != np.floor(p):
            p = np.floor(p)
 
        if np.iscomplex(p) or p < 1:
            raise ValueError(''p must be at least 1 and real'')
 
        pn = np.sum(np.absolute(self.x) ** p) ** (1. / p)
 
        return pn

def __str__(self):
        # create local dictionary to convert node-numbers to node-labels
        num2node_label = {num: name for name, num in self.node_label2num.items()}
 
        # build message to print
        msg = ''------------------------\\n''
        msg += ''    SpicePy.Network:\\n''
        msg += ''------------------------\\n''
        for ele, nodes, val in zip(self.names, self.nodes, self.values):
            # if val is a list --> ele is a transient source
            if isinstance(val, list):
                if self.source_type[ele] == ''pwl'':
                    fmt = "{} {} {} {}(" + "{} " * (len(val[0]) - 1) + "{})\\n"
                    msg += fmt.format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], self.source_type[ele], *val[0])
                else:
                    fmt = "{} {} {} {}(" + "{} " * (len(val) - 1) + "{})\\n"
                    msg += fmt.format(ele, *val)
            # if val is complex --> ele is a phasor
            elif np.iscomplex(val):
                msg += "{} {} {} {} {}\\n".format(ele, np.abs(val), np.angle(val) * 180/np.pi)
            # if ele is C or L
            elif ele[0].upper() == ''C'' or ele[0].upper() == ''L'':
                # check if an i.c. is present and print it
                if ele in self.IC:
                    msg += "{} {} {} {} ic={}\\n".format(ele, val, self.IC[ele])
                # otherwise...
                else:
                    msg += "{} {} {} {}\\n".format(ele, val)
            # otherwise...general case -->  ele n+ n- val
            else:
                msg += "{} {} {} {}\\n".format(ele, val)
 
        # add analysis
        msg += " ".join(self.analysis) + ''\\n''
 
        # if a plot command is present,add it
        if self.plot_cmd is not None:
            msg += self.plot_cmd + ''\\n''
 
        # add number of nodes (reference node is included) and number of branches
        msg += ''------------------------\\n''
        msg += ''* number of nodes {}\\n''.format(self.node_num + 1)
        msg += ''* number of branches {}\\n''.format(len(self.names))
        msg += ''------------------------\\n''
 
        return msg

def find_zero_crossing_quadratic(x1, y1, x2, y2, x3, y3, eps = 1.0):
        """ Find zero crossing using quadratic approximation along 1d line"""
        # compute coords along 1d line
        v = x2 - x1
        v = v / np.linalg.norm(v)
        if v[v!=0].shape[0] == 0:
            logging.error(''Difference is 0. Probably a bug'')
 
        t1 = 0
        t2 = (x2 - x1)[v!=0] / v[v!=0]
        t2 = t2[0]
        t3 = (x3 - x1)[v!=0] / v[v!=0]
        t3 = t3[0]
 
        # solve for quad approx
        x1_row = np.array([t1**2, t1, 1])
        x2_row = np.array([t2**2, t2, 1])
        x3_row = np.array([t3**2, t3, 1])
        X = np.array([x1_row, x2_row, x3_row])
        y_vec = np.array([y1, y3])
        try:
            w = np.linalg.solve(X, y_vec)
        except np.linalg.LinAlgError:
            logging.error(''Singular matrix. Probably a bug'')
            return None
 
        # get positive roots
        possible_t = np.roots(w)
        t_zc = None
        for i in range(possible_t.shape[0]):
            if possible_t[i] >= 0 and possible_t[i] <= 10 and not np.iscomplex(possible_t[i]):
                t_zc = possible_t[i]
 
        # if no positive roots find min
        if np.abs(w[0]) < 1e-10:
            return None
 
        if t_zc is None:
            t_zc = -w[1] / (2 * w[0])
 
        if t_zc < -eps or t_zc > eps:
            return None
 
        x_zc = x1 + t_zc * v
        return x_zc

def run(X,y1,y2,dmax=None):
 
    D,N = X.shape
    y1_unique,J1 = np.unique(y1, return_inverse=True)
    ny1 = y1_unique.size
 
    y2_unique,J2 = np.unique(y2, return_inverse=True)
    ny2 = y2_unique.size
 
    Y1 = np.zeros((ny1,N))
    Y2 = np.zeros((ny2,N))
    Y1[J1,range(N)] = 1.
    Y2[J2,range(N)] = 1.
 
    XY2 = np.dot(X,Y2.T) # D x ny2
    XY2Y2X = np.dot(XY2,XY2.T) # D x D
    XX = np.dot(X,X.T) # D x D
 
    P = np.zeros((D,0))
    Sj = np.dot(X,Y1.T) #  D x ny1
 
    if dmax==None:
        dmax = D
 
    for d in range(dmax):
        if d>0: 
            invPP = np.linalg.pinv(np.dot(P.T,P))
            Sj -= np.dot(np.dot(np.dot(P,invPP),P.T),Sj) 
 
        C = np.dot(Sj,Sj.T) - XY2Y2X
        C = 0.5*(C+C.T)
        dd,E = scipy.linalg.eigh(C,eigvals=(D-1,D-1)) # ascending order
 
        assert np.isnan(dd).any()==False
        assert np.iscomplex(dd).any()==False
        #dd = dd[::-1] #
        #E = E[:,::-1]
        wj = E#E[:,0] # D x 1
        pj = np.dot(XX,wj) / np.dot(np.dot(wj.T,XX),wj) # D x 1
        P = np.hstack((P,pj.reshape((D,1)))) #  D x d
 
 
    #P = P/np.tile(np.sqrt((P**2).sum(axis=0,keepdims=True)),(D,1))
    #% They need not be orthogonal.
    return P

def run(X,y2):
    # function [E,dd] = privacyLDA(X,y1,y2)
    # %max_W0 tr(W0''*C1*W0) - tr(W0''*C2*W0)
    D,N = X.shape
    y1_unique = np.unique(y1)
    #print y1_unique
    #[y1_unique,~,J1] = unique(y1);
    #ny1 = y1_unique.size
 
    y2_unique = np.unique(y2)
    #print y2_unique
    #[y2_unique,J2] = unique(y2);
    #ny2 = y2_unique.size
 
    C1 = np.zeros((D,D))
    C2 = np.zeros((D,D))
    mu = X.mean(axis=1).reshape((D,1))
    #print mu.shape
 
    for k in np.nditer(y1_unique):
        indk = np.where(y1==k)[0]        
        muk = X[:,indk].mean(axis=1).reshape((D,1))
        #muk -= np.kron(np.ones((1,len(indk))),mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)'';
        C1 = C1 + len(indk)*np.dot(muk-mu,(muk-mu).T)
 
    for k in np.nditer(y2_unique):
        indk = np.where(y2==k)[0]        
        muk = X[:,mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)'';
        C2 = C2 + len(indk)*np.dot(muk-mu,(muk-mu).T)
 
    C1 = C1 + 1e-8*np.trace(C1)*np.eye(D)# 
    C2 = C2 + 1e-8*np.trace(C2)*np.eye(D)#
    C1 = 0.5*(C1+C1.T)#;% + 1E-8*trace(C1)*eye(D); 
    C2 = 0.5*(C2+C2.T)#;% + 1E-8*trace(C2)*eye(D);
 
 
    dd,E = scipy.linalg.eigh(C1,C2) # ascending order
 
    #print dd.shape
    #print E.shape
 
    assert np.isnan(dd).any()==False
    assert np.iscomplex(dd).any()==False
    #[dd,ind] = sort(diag(dd),''descend''); 
    #print dd
    dd = dd[::-1] #
    E = E[:,::-1]
 
    E = E/np.tile(np.sqrt((E**2).sum(axis=0,keepdims=True)),(D,1))
    #% They need not be orthogonal.
 
    #print dd.shape
    #print E.shape
 
    return (E,dd)

def cwplot(fb_est,rx,t,fs:int,fn) -> None:
#%% time
    fg,axs = subplots(1,2,figsize=(12,6))
    ax = axs[0]
    ax.plot(t, rx.T.real)
    ax.set_xlabel(''time [sec]'')
    ax.set_ylabel(''amplitude'')
    ax.set_title(''Noisy,jammed receive signal'')
    #%% periodogram
    if DTPG >= (t[-1]-t[0]):
        dt = (t[-1]-t[0])/4
    else:
        dt = DTPG
 
    dtw = 2*dt #  seconds to window
    tstep = ceil(dt*fs)
    wind = ceil(dtw*fs);
    nfft = zeropadfactor*wind
 
    f,Sraw = signal.welch(rx.ravel(), fs,
                          nperseg=wind,noverlap=tstep,nfft=nfft,
                          return_onesided=False)
 
    if np.iscomplex(rx).any():
        f = np.fft.fftshift(f); Sraw = np.fft.fftshift(Sraw)
 
    ax=axs[1]
    ax.plot(f,Sraw,''r'',label=''raw signal'')
 
    fc_est = f[Sraw.argmax()]
 
    #ax.set_yscale(''log'')
    ax.set_xlim([fc_est-200,fc_est+200])
    ax.set_xlabel(''frequency [Hz]'')
    ax.set_ylabel(''amplitude'')
    ax.legend()
 
    esttxt=''''
 
    if fn is None: # simulation
        ax.axvline(ft+fb0,color=''red'',linestyle=''--'',label=''true freq.'')
        esttxt += f''true: {ft+fb0} Hz ''
 
    for e in fb_est:
        ax.axvline(e,color=''blue'',label=''est. freq.'')
 
    esttxt += '' est: '' + str(fb_est) +'' Hz''
 
    ax.set_title(esttxt)

def cw_est(rx, fs:int, Ntone:int, method:str=''esprit'', usepython=False, useall=False):
    """
    estimate beat frequency using subspace frequency estimation techniques.
    This is much faster in Fortran,but to start using Python alone doesn''t require compiling Fortran.
 
    ESPRIT and RootMUSIC are two popular subspace techniques.
 
    Matlab''s rootmusic is a far inferior FFT-based method with very poor accuracy vs. my implementation.
    """
    assert isinstance(method,str)
    method = method.lower()
 
    tic = time()
    if method == ''esprit'':
#%% ESPRIT
        if rx.ndim == 2:
            assert usepython,''Fortran not yet configured for multi-pulse case''
            Ntone *= 2
 
 
        if usepython or (Sc is None and Sr is None):
            print(''Python ESPRIT'')
            fb_est, sigma = esprit(rx, Ntone, Nblockest, fs)
        elif np.iscomplex(rx).any():
            print(''Fortran complex64 ESPRIT'')
            fb_est, sigma = Sc.subspace.esprit(rx,Ntone,Nblockest,fs)
        else: # real signal
            print(''Fortran float32 ESPRIT'')
            fb_est, sigma = Sr.subspace.esprit(rx,fs)
 
        fb_est = abs(fb_est)
#%% ROOTMUSIC
    elif method == ''rootmusic'':
        fb_est, sigma = rootmusic(rx,fs)
    else:
        raise ValueError(f''unkNown estimation method: {method}'')
    print(f''computed via {method} in {time()-tic:.1f} seconds.'')
#%% improvised process for CW only without notch filter
    # assumes first two results have largest singular values (from SVD)
    if not useall:
        i = sigma > 0.001 # arbitrary
        fb_est = fb_est[i]
        sigma  = sigma[i]
 
#        if fb_est.size>1:
#            ii = np.argpartition(sigma,Ntone-1)[:Ntone-1]
#            fb_est = fb_est[ii]
#            sigma = sigma[ii]
 
 
    return fb_est, sigma