def minZerrKernSHG_naive(self):
        Et0 = self.Et_cla.get()
        Esig = self.Esig_t_tau_p_cla.get()
        dZ = self.dZ_cla.get()
        N = Esig.shape[0]

        mx = 0.0
        X = np.zeros(5)

        for t in range(N):
            for tau in range(N):
                T = np.abs(Esig[tau, t]) ** 2
                if mx < T:
                    mx = T
                tp = t - (tau - N / 2)
                if tp >= 0 and tp < N:
                    dZdZ = dZ[t] * dZ[tp]
                    dZE = dZ[t] * Et0[tp] + dZ[tp] * Et0[t]
                    DEsig = Et0[t] * Et0[tp] - Esig[tau, t]

                    X[0] += np.abs(dZdZ) ** 2
                    X[1] += 2.0 * np.real(dZE * np.conj(dZdZ))
                    X[2] += 2.0 * np.real(DEsig * np.conj(dZdZ)) + np.abs(dZE) ** 2
                    X[3] += 2.0 * np.real(DEsig * np.conj(dZE))
                    X[4] += np.abs(DEsig) ** 2
        T = N * N * mx
        X[0] = X[0] / T
        X[1] = X[1] / T
        X[2] = X[2] / T
        X[3] = X[3] / T
        X[4] = X[4] / T

        root.debug("".join(("Esig_t_tau_p norm max: ", str(mx))))

        return X

def Dmat(alpha,xAve,pAve):

    D = np.zeros((Nb,Nb),dtype=complex)

    V0, V1, ak  = Hessian(xAve) 

    # time derivative 
    dq = pAve / am
    dp = - V1
    da = (-1.0j/am * alpha**2 + 1.j * ak) 
    
    a = alpha.real 

    for k in range(Nb):
        D[k,k] = - 1j*da.imag/2./a * (float(k) + 0.5) - 1j * pAve**2/am    
    
    for k in range(1,Nb):
        D[k-1,k] = np.sqrt(float(k)/2./a) * ( - np.conj(alpha) * dq + 1j * dp) 
        D[k,k-1] = np.sqrt(float(k)/2./a) * ( alpha * dq + 1j * dp )
    
    if Nb > 2: 
        for k in range(2,Nb):
            D[k-2,k] = np.conj(da)/2./a * np.sqrt(float(k*(k-1)))/2.0 
            D[k,k-2] = - da/2./a * np.sqrt(float(k*(k-1)))/2.0
            


    #D[0,0] = - 1j * pAve**2 / am 
    #D[1,1] = - 1j * pAve**2 / am 
    #D[0,1] = - (np.conj(alpha)*pAve/am + 1j * V1) / np.sqrt(2.*alpha.real)
    #D[1,0] = (alpha * pAve / am - 1j * V1) / np.sqrt(2.*alpha.real)

    return D 

def childGeneration(N, particles, tree, i):
    maxInCell = 3
    p, t = noOfParticlesInside(particles, N[i])
    N[i].particleIndex = t
    if p > maxInCell:
        t = len(N)
        tree[i] = [t, t + 1, t + 2, t + 3]
        c = N[i].center
        v = N[i].vertex
        N.append(Node((c - (v - c) / 2), c, t))
        N.append(Node((c + (v - c) / 2), v, t + 1))
        N.append(Node(c + (conj(v - c) / 2), real(v) + 1j * imag(c), t + 2))
        N.append(Node(c - (conj(v - c) / 2), real(c) + 1j * imag(v), t + 3))
        tree.append([])
        tree.append([])
        tree.append([])
        tree.append([])
        for ch in tree[i]:
            childGeneration(N, particles, tree, ch)

        N[i].aj = zeros(noOfTerms) * 1j
        for j in range(1, noOfTerms + 1):
            for k in range(1, j + 1):
                for chi in tree[i]:
                    N[i].aj[j - 1] += (
                        N[chi].aj[k - 1] * combitorial(j - 1, k - 1) * pow((-N[i].center + N[chi].center), j - k)
                    )

    else:
        N[i].aj = zeros(noOfTerms) * 1j
        for j in range(1, noOfTerms + 1):
            for tempIdx in t:
                N[i].aj[j - 1] += particles[tempIdx].strength * pow((particles[tempIdx].xy - N[i].center), (j - 1))

    def mfunc(uv, p, d, f):
        global curtime
        crd,t,(i,j) = p
        bl = a.miriad.ij2bl(i,j)
        if bl in conj_bls: d = n.conj(d)
        if is_run1(t):
            #if i == 2 or j == 2: return p, None, None
            d = n.conj(d)
            i,j = rewire_run1[i], rewire_run1[j]
            uvo['pol'] = a.miriad.str2pol['xx']
        elif is_run2(t):
            #if i == 2 or j == 2: return p, None, None
	    d = n.conj(d)
	    i,j = rewire_run2[i], rewire_run2[j]
            uvo['pol'] = a.miriad.str2pol['xx']
        else: return p, None, None
        if i > j: i,j,d = j,i,n.conj(d)
        if curtime != t:
            #if is_run1(t): print 'Processing as run 1'
            #elif is_run2(t): print 'Processing as run 2'
            curtime = t
            aa.set_jultime(t)
            uvo['lst'] = aa.sidereal_time()
            uvo['ra'] = aa.sidereal_time()
            uvo['obsra'] = aa.sidereal_time()
        p = crd,t,(i,j)
        return p,d,f

def get_modal_volume(sim, box=None, dft_cell=None, nf=0):

    Exyz=[mp.Ex, mp.Ey, mp.Ez]

    if dft_cell is None:
      (Ex,Ey,Ez) = [sim.get_array(vol=box, component=c, cmplx=True) for c in Exyz]
      (X,Y,Z,W)  = sim.get_array_metadata(vol=box)
      Eps        = sim.get_array(vol=box, component=mp.Dielectric)
    else:
      (Ex,Ey,Ez) = [sim.get_dft_array(dft_cell, c, nf) for c in Exyz]
      (X,Y,Z,W)  = sim.get_dft_array_metadata(dft_cell=dft_cell)
      # slightly annoying: we need an epsilon array with empty dimensions collapsed,
      #  something not currently provided by any C++ function; for now just
      #  create it via brute-force python loop, although this could get slow.
      Eps=np.zeros(0)
      for x in X:
        for y in Y:
          for z in Z:
            Eps=np.append(Eps,sim.get_epsilon_point(mp.Vector3(x,y,z)))
      Eps=np.reshape(Eps,np.shape(W))

    EpsE2 = np.real(Eps*(np.conj(Ex)*Ex + np.conj(Ey)*Ey + np.conj(Ez)*Ez))
    num   = np.sum(W*EpsE2)
    denom = np.max(EpsE2)

    return num/denom if denom!=0.0 else 0.0

def visualize_dft_flux(sim, superpose=True, flux_cells=[],
                       options=None, nf=0):

    if not mp.am_master():
        return
    options=options if options else def_flux_options

    # first pass to get arrays of poynting flux strength for all cells
    if len(flux_cells)==0:
        flux_cells=[cell for cell in sim.dft_objects if is_flux_cell(cell)]
    flux_arrays=[]
    for cell in flux_cells:    # first pass to compute flux data
        (x,y,z,w,c,EH)=unpack_dft_cell(sim,cell,nf=nf)
        flux_arrays.append( 0.25*np.real(w*(np.conj(EH[0])*EH[3] - np.conj(EH[1])*EH[2])) )

    # second pass to plot
    for n, cell in enumerate(flux_cells):    # second pass to plot
        if superpose==False:
            if n==0:
                plt.figure()
                plt.title('Poynting flux')
            plt.subplot(len(flux_cells),1,n)
            plt.gca().set_title('Flux cell {}'.format(n))
        cn,sz=mp.get_center_and_size(cell.where)
        max_flux=np.amax([np.amax(fa) for fa in flux_arrays])
        plot_data_curves(sim, center=cn, size=sz, data=[flux_arrays[n]],
                         superpose=superpose, options=options,
                         labels=['flux through cell {}'.format(n)],
                         dmin=-max_flux,dmax=max_flux)

def alignShells(alm, blm, lmax, Lgrid):

    L = lmax+1
    euler = sampling_grid(Lgrid)
    N = np.size(euler, 0)

    Corr = np.zeros(N, dtype=complex)
    for l in range(L):
        for m in range(-l,l+1):
            for n in range(-l,l+1):
                dlmn = wignerD(l,m,n)(euler[:,0],euler[:,1],euler[:,2])

                a = hp2lm(alm, l, m, lmax)
                b = hp2lm(blm,l,n,lmax)
                Corr+= a*np.conj(b)*np.conj(dlmn)

    x = np.where(Corr==max(Corr))
    x = x[0]

    # now we rotate back blm         
    blmr = rotatealm(blm,lmax,euler[x,0][0],euler[x,1][0],euler[x,2][0])

    hpblmr = stand2hp(blmr,lmax)

    return hpblmr

 def bartlet(self, sv): 
   ''' Bartlet's estimator for DOA. Use `argmax`. ''' 
   V = self.signal_vector 
   G = sv.steering_vectors[self.site_id] 
   self.bearing = sv.bearings[self.site_id]
   left_half = np.dot(V, np.conj(np.transpose(G))) 
   return np.real(left_half * np.conj(left_half)) 

    def test_adjoints(exp, freqs):

        deltas = exp.m0.mesh.deltas

        m1_ = exp.m1.asarray()

        lindata = exp.lin_results["simdata"]
        dhat = exp.dhat

        adjmodel = exp.adj_results["imaging_condition"].asarray()

        temp_data_prod = 0.0
        for nu in freqs:
            temp_data_prod += np.dot(lindata[nu].reshape(dhat[nu].shape), np.conj(dhat[nu]))

        pt1 = temp_data_prod

        pt2 = np.dot(m1_.T, np.conj(adjmodel)).squeeze() * np.prod(deltas)

        print "{0}: ".format(exp.name)
        print "<Fm1, d>             = {0: .4e} ({1:.4e})".format(pt1, np.linalg.norm(pt1))
        print "<m1, F*d>            = {0: .4e} ({1:.4e})".format(pt2, np.linalg.norm(pt2))
        print "<Fm1, d> - <m1, F*d> = {0: .4e} ({1:.4e})".format(pt1 - pt2, np.linalg.norm(pt1 - pt2))

        print "Relative error       = {0: .4e}\n".format(np.linalg.norm(pt1 - pt2) / np.linalg.norm(pt1))

def test1():
    print("*** MPS tests started ***")
    (N,chi,d) = (7,10,2)
    A = randomMPS(N,chi,d)
    state = getState(A)
    state = state/np.sqrt(np.dot(np.conj(state),state))
    prod = np.dot(np.conj(state),state)
    approxA = getMPS(state,2)
    approxState = getState(approxA)
    approxProd = np.dot(np.conj(approxState),approxState)
    relErr = approxProd/prod - 1
    S = entropy(state)
    print("State total %d elements"%state.size)
    print("MPS total %d elements"%A.size)
    print("(N,chi,d) = (%d,%d,%d)"%(N,chi,d))
    print("Expected:        (%f,%f)"%polar(prod))
    print("SVD:             (%f,%f)"%polar(innerProduct(approxA,approxA)))
    print("Product:         (%f,%f)"%polar(approxProd))
    print("Relative error:  %f"%np.absolute(relErr))
    print("Entropy:         %f"%S)
    print("")
    # state = np.ones(d**N)/np.sqrt(2)**N
    # state = np.zeros(2**10)
    # state[0] = 1/np.sqrt(2)
    # state[-1] = 1/np.sqrt(2)
    state = np.random.rand(d**N)
    state = state/np.linalg.norm(state)
    
    mps = getMPS(state,4)
    print("Expected: (%f,%f)"%polar(np.inner(state,state)))
    print("MPS:      (%f,%f)"%polar(innerProduct(mps,mps)))
    print("*** MPS tests finished ***\n")

def par_xstep(i):
    r"""Minimise Augmented Lagrangian with respect to
    :math:`\mathbf{x}_{G_i}`, one of the disjoint problems of optimizing
    :math:`\mathbf{x}`.

    Parameters
    ----------
    i : int
      Index of grouping to update

    """
    global mp_X
    global mp_DX
    YU0f = sl.rfftn(mp_Y0[[i]] - mp_U0[[i]], mp_Nv, mp_axisN)
    YU1f = sl.rfftn(mp_Y1[mp_grp[i]:mp_grp[i+1]] -
                    1/mp_alpha*mp_U1[mp_grp[i]:mp_grp[i+1]], mp_Nv, mp_axisN)
    if mp_Cd == 1:
        b = np.conj(mp_Df[mp_grp[i]:mp_grp[i+1]]) * YU0f + mp_alpha**2*YU1f
        Xf = sl.solvedbi_sm(mp_Df[mp_grp[i]:mp_grp[i+1]], mp_alpha**2, b,
                            mp_cache[i], axis=mp_axisM)
    else:
        b = sl.inner(np.conj(mp_Df[mp_grp[i]:mp_grp[i+1]]), YU0f,
                     axis=mp_C) + mp_alpha**2*YU1f
        Xf = sl.solvemdbi_ism(mp_Df[mp_grp[i]:mp_grp[i+1]], mp_alpha**2, b,
                              mp_axisM, mp_axisC)
    mp_X[mp_grp[i]:mp_grp[i+1]] = sl.irfftn(Xf, mp_Nv,
                                            mp_axisN)
    mp_DX[i] = sl.irfftn(sl.inner(mp_Df[mp_grp[i]:mp_grp[i+1]], Xf,
                                  mp_axisM), mp_Nv, mp_axisN)

def test_solve_fredholm_reconstr_ac():
    """
        here we see that the reconstruction quality is independent of the integration weights

        differences occur when checking validity of the interpolated time continuous Fredholm equation
    """
    _WC_ = 2
    def lac(t):
        return np.exp(- np.abs(t) - 1j*_WC_*t)
    t_max = 10
    tol = 2e-10
    for ng in range(11,500,30):
        t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng)
        r = lac(t.reshape(-1,1)-t.reshape(1,-1))
        _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w)
        _eig_vec_ast = np.conj(_eig_vec)  # (N_gp, N_ev)
        tmp = _eig_val.reshape(1, -1) * _eig_vec  # (N_gp, N_ev)
        recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1]))
        rd = np.max(np.abs(recs_bcf - r) / np.abs(r))
        assert rd < tol, "rd={} >= {}".format(rd, tol)

        t, w = sp.method_kle.get_simpson_weights_times(t_max, ng)
        r = lac(t.reshape(-1, 1) - t.reshape(1, -1))
        _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w)
        _eig_vec_ast = np.conj(_eig_vec)  # (N_gp, N_ev)
        tmp = _eig_val.reshape(1, -1) * _eig_vec  # (N_gp, N_ev)
        recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1]))
        rd = np.max(np.abs(recs_bcf - r) / np.abs(r))
        assert rd < tol, "rd={} >= {}".format(rd, tol)

def calc_stokes(Ex,Ey,s=None):
    
    if len(Ex) != len(Ey):
        raise ValueError('Ex and Ey dimentions do not match')
        
    if s is None:
        s = np.arange(len(Ex))
    
    Ex_ = np.conj(Ex)
    Ey_ = np.conj(Ey)
    
    Jxx = Ex * Ex_
    Jxy = Ex * Ey_
    Jyx = Ey * Ex_
    Jyy = Ey * Ey_
    
    del (Ex_,Ey_)
    
    S = StokesParameters()
    S.sc = s
    S.s0 = real(Jxx + Jyy)
    S.s1 = real(Jxx - Jyy)
    S.s2 = real(Jxy + Jyx)
    S.s3 = real(1j * (Jyx - Jxy))
    
    return S

def deconstruct_single_qubit_matrix_into_angles(
        mat: np.ndarray) -> Tuple[float, float, float]:
    """Breaks down a 2x2 unitary into more useful ZYZ angle parameters.

    Args:
        mat: The 2x2 unitary matrix to break down.

    Returns:
        A tuple containing the amount to phase around Z, then rotate around Y,
        then phase around Z (all in radians).
    """
    # Anti-cancel left-vs-right phase along top row.
    right_phase = cmath.phase(mat[0, 1] * np.conj(mat[0, 0])) + math.pi
    mat = np.dot(mat, _phase_matrix(-right_phase))

    # Cancel top-vs-bottom phase along left column.
    bottom_phase = cmath.phase(mat[1, 0] * np.conj(mat[0, 0]))
    mat = np.dot(_phase_matrix(-bottom_phase), mat)

    # Lined up for a rotation. Clear the off-diagonal cells with one.
    rotation = math.atan2(abs(mat[1, 0]), abs(mat[0, 0]))
    mat = np.dot(_rotation_matrix(-rotation), mat)

    # Cancel top-left-vs-bottom-right phase.
    diagonal_phase = cmath.phase(mat[1, 1] * np.conj(mat[0, 0]))

    # Note: Ignoring global phase.
    return right_phase + diagonal_phase, rotation * 2, bottom_phase

    def mouse_move(self, event):
        """ Handle mouse movement, redraw pzmap while drag/dropping """
        if (self.move_zero != None and
            event.xdata != None and 
            event.ydata != None):

            if (self.index1 == self.index2):
                # Real pole/zero
                new = event.xdata
                if (self.move_zero == True):
                    self.zeros[self.index1] = new
                elif (self.move_zero == False):
                    self.poles[self.index1] = new
            else:
                # Complex poles/zeros
                new = event.xdata + 1.0j*event.ydata
                if (self.move_zero == True):
                    self.zeros[self.index1] = new
                    self.zeros[self.index2] = conj(new)
                elif (self.move_zero == False):
                    self.poles[self.index1] = new
                    self.poles[self.index2] = conj(new)
            tfcn = None
            if (self.move_zero == True):
                self.tfi.set_zeros(self.zeros)
                tfcn = self.tfi.get_tf()
            elif (self.move_zero == False):
                self.tfi.set_poles(self.poles)
                tfcn = self.tfi.get_tf()
            if (tfcn != None):
                self.draw_pz(tfcn)
                self.canvas_pzmap.draw()

def autocorr_fft(signal, axis = -1):
    """Return full autocorrelation along specified axis. Use fft
    for computation."""
    if N.ndim(signal) == 0:
        return signal
    elif signal.ndim == 1:
        n       = signal.shape[0]
        nfft    = int(2 ** nextpow2(2 * n - 1))
        lag     = n - 1
        a       = fft(signal, n = nfft, axis = -1)
        au      = ifft(a * N.conj(a), n = nfft, axis = -1)
        return N.require(N.concatenate((au[-lag:], au[:lag+1])), dtype = signal.dtype)
    elif signal.ndim == 2:
        n       = signal.shape[axis]
        lag     = n - 1
        nfft    = int(2 ** nextpow2(2 * n - 1))
        a       = fft(signal, n = nfft, axis = axis)
        au      = ifft(a * N.conj(a), n = nfft, axis = axis)
        if axis == 0:
            return N.require(N.concatenate( (au[-lag:], au[:lag+1]), axis = axis), \
                    dtype = signal.dtype)
        else:
            return N.require(N.concatenate( (au[:, -lag:], au[:, :lag+1]), 
                        axis = axis), dtype = signal.dtype)
    else:
        raise RuntimeError("rank >2 not supported yet")

def deconvolve(signal, HRFs, SNRs):
	"""Deconvolve signal using Wiener deconvolution."""
	H = fft(HRFs, len(signal), axis=0)
	wiener_filter = np.conj(H) / (H * np.conj(H) + 1 / SNRs**2)
	deconvolved = np.real(ifft(wiener_filter * fft(signal, axis=0), axis=0))

	return deconvolved

def levup(acur, knxt, ecur=None):
    """LEVUP  One step forward Levinson recursion

    :param acur:
    :param knxt:
    :return: 
        * anxt the P+1'th order prediction polynomial based on the P'th order prediction polynomial, acur, and the 
          P+1'th order reflection coefficient, Knxt.
        * enxt the P+1'th order prediction  prediction error, based on the P'th order prediction error, ecur.  


    :References:  P. Stoica R. Moses, Introduction to Spectral Analysis  Prentice Hall, N.J., 1997, Chapter 3.
    """
    if acur[0] != 1:
        raise ValueError("At least one of the reflection coefficients is equal to one.")
    acur = acur[1:]  #  Drop the leading 1, it is not needed

    # Matrix formulation from Stoica is used to avoid looping
    anxt = numpy.concatenate((acur, [0])) + knxt * numpy.concatenate((numpy.conj(acur[-1::-1]), [1]))

    enxt = None
    if ecur != None:
        # matlab version enxt = (1-knxt'.*knxt)*ecur
        enxt = (1.0 - numpy.dot(numpy.conj(knxt), knxt)) * ecur

    anxt = numpy.insert(anxt, 0, 1)

    return anxt, enxt

def vis_cal(visdata,vismodel,max_iter=2000,threshold=0.0001):

	nchan = visdata.shape[-1]
	nant = visdata.shape[0]
	gains=NP.ones((nant,nchan),dtype=NP.complex64)

	# set up bookkeeping
	ant1 = NP.arange(nant)

	chi_history = NP.zeros((max_iter,nchan))

	tries = 0.0
	change = 100.0
	A=NP.zeros((nant**2,nant),dtype=NP.complex64) # matrix for minimization
	ind1 = NP.arange(nant**2)
	ind2 = NP.repeat(NP.arange(nant),nant)
	ind3 = NP.tile(NP.arange(nant),nant)
	for fi in xrange(nchan):
		tempgains = gains[:,fi].copy()
		tempdata = visdata[:,:,fi].reshape(-1)
		tempmodel = vismodel[:,:,fi].reshape(-1)
		while (tries < max_iter) and (change > threshold):
			chi_history[tries,fi] = NP.sum(tempdata-NP.outer(tempgains,NP.conj(tempgains))*tempmodel)
			prevgains = tempgains.copy()
			A[ind1,ind2] = tempmodel*NP.conj(prevgains[ind3])
			tempgains = NP.linalg.lstsq(A, tempdata)[0]
			change = NP.median(NP.abs(tempgains-prevgains)/NP.abs(prevgains))
			tries += 1
		if tries == max_iter:
			print 'Warning! Vis calibration failed to converge. Continuing'
		gains[:,fi] = tempgains.copy()

	return gains

 def right_sweep(self,R,M):
     """
     :param R: The right operator
     :param M:
     :return:
     """
     H=np.tensordot(self.MPO[0],R[1],axes=(3,1))
     H=np.squeeze(H)
     H=np.transpose(H,[0,2,1,3])
     d0,d1,d2,d3=H.shape
     H=np.reshape(H,[d0*d1,d2*d3])
     w,v=ssl.eigsh(H,which='SA',k=1,maxiter=5000)
     v=np.reshape(v,[self.d,1,d0*d1//self.d])
     l,u=self.left_cannonical(v)
     M[0]=l
     L=[[] for i in range(len(R))]
     L[0]=np.tensordot(l,self.MPO[0],axes=(0,0))
     L[0]=np.tensordot(L[0],np.conj(l),axes=(2,0))
     L[0]=np.transpose(L[0],[0,2,4,1,3,5])
     for i in range(1,len(R)-1):
         H=np.tensordot(self.MPO[i],R[i+1],axes=(3,1))
         H=np.tensordot(L[i-1],H,axes=(4,2))
         H=H.squeeze()
         H=np.transpose(H,[2,0,4,3,1,5])
         d1,d2,d3,d4,d5,d6=H.shape
         H=np.reshape(H,[d1*d2*d3,d1*d2*d3])
         w,v=ssl.eigsh(H,which='SA',k=1,maxiter=5000)
         v = np.reshape(v, [d1, d2, d3])
         l, u = self.left_cannonical(v)
         M[i] = l
         Li = np.tensordot(L[i-1],l,axes=(3,1))
         Li = np.tensordot(Li, self.MPO[i],axes=([5,3],[0,2]))
         L[i] = np.tensordot(Li, np.conj(l),axes=([3,5],[1,0]))
     M[-1]=np.tensordot(u,M[-1],axes=(1,1))
     return L,M