def is_finite(x):
    try:
        if type(x) == str:
            return False
        return np.array([np.isreal(e) and np.isfinite(e) for e in x], dtype=bool)
    except TypeError:
        return np.isreal(x) and np.isfinite(x)

def is_nan(x):
    try:
        if type(x) == str:
            raise TypeError()
        return np.array([e == NA_character_ or (np.isreal(e) and np.isnan(e)) for e in x], dtype=bool)
    except TypeError:
        return x == NA_character_ or (np.isreal(x) and np.isnan(x))

def test_multib0_dsi():
    data, gtab = dsi_voxels()
    # Create a new data-set with a b0 measurement:
    new_data = np.concatenate([data, data[..., 0, None]], -1)
    new_bvecs = np.concatenate([gtab.bvecs, np.zeros((1, 3))])
    new_bvals = np.concatenate([gtab.bvals, [0]])
    new_gtab = gradient_table(new_bvals, new_bvecs)
    ds = DiffusionSpectrumModel(new_gtab)
    sphere = get_sphere('repulsion724')
    dsfit = ds.fit(new_data)
    pdf = dsfit.pdf()
    dsfit.odf(sphere)
    assert_equal(new_data.shape[:-1] + (17, 17, 17), pdf.shape)
    assert_equal(np.alltrue(np.isreal(pdf)), True)

    # And again, with one more b0 measurement (two in total):
    new_data = np.concatenate([data, data[..., 0, None]], -1)
    new_bvecs = np.concatenate([gtab.bvecs, np.zeros((1, 3))])
    new_bvals = np.concatenate([gtab.bvals, [0]])
    new_gtab = gradient_table(new_bvals, new_bvecs)
    ds = DiffusionSpectrumModel(new_gtab)
    dsfit = ds.fit(new_data)
    pdf = dsfit.pdf()
    dsfit.odf(sphere)
    assert_equal(new_data.shape[:-1] + (17, 17, 17), pdf.shape)
    assert_equal(np.alltrue(np.isreal(pdf)), True)

def findSigT(xIn, yIn):
    ##
    ## Added to check their sig test. Yes I know it is paranoid rjel
    ##   regPoints are the number of points going into the
    ##    regression calculation.
    ## Check that there are some points not masked out first. Else you
    ##        get horrible errors.
    if np.any(np.isreal(((xIn <= 1.0) & (np.isreal(yIn))))):
        rPoints = int(((xIn <= 1.0) & np.isreal(yIn)).sum())
    else:
        rPoints = 0
    ##    print( 'There are ', regPoints,' points in the calculation')
    ##
    ## Find the value in the ttest table corresponding to n and
    ##    the p value of interest. The table I downloaded had 1-100
    ##        points then more than 100
    if (rPoints < 103) & (rPoints > 2):
        nIndex = rPoints - 2
    ##
    ## if there are less than 2 points don't bother
    elif rPoints <= 2:
        return 3000.0, rPoints
    else:
        nIndex = 101
    ## The first line of the t-table file has the p values in it
    ##  so we can search a split first line for the column
    ##   that contains a given p-value
    ##
    pIndex = tTable[0].index(str(1.0 - p.pVal))
    sigLocal = tTable[nIndex][pIndex]

    return sigLocal, rPoints

def filter_params(io, rsh, rs, ee, isc):
    # Function filter_params identifies bad parameter sets. A bad set contains
    # Nan, non-positive or imaginary values for parameters; Rs > Rsh; or data
    # where effective irradiance Ee differs by more than 5% from a linear fit
    # to Isc vs. Ee

    badrsh = np.logical_or(rsh < 0., np.isnan(rsh))
    negrs = rs < 0.
    badrs = np.logical_or(rs > rsh, np.isnan(rs))
    imagrs = ~(np.isreal(rs))
    badio = np.logical_or(~(np.isreal(rs)), io <= 0)
    goodr = np.logical_and(~badrsh, ~imagrs)
    goodr = np.logical_and(goodr, ~negrs)
    goodr = np.logical_and(goodr, ~badrs)
    goodr = np.logical_and(goodr, ~badio)

    matrix = np.vstack((ee / 1000., np.zeros(len(ee)))).T
    eff = np.linalg.lstsq(matrix, isc)[0][0]
    pisc = eff * ee / 1000
    pisc_error = np.abs(pisc - isc) / isc
    # check for departure from linear relation between Isc and Ee
    badiph = pisc_error > .05

    u = np.logical_and(goodr, ~badiph)
    return u

def bezier_curve_y_out(shift_angle, P0, P1, P2, P3, second_of_day=None):
    '''
    For a cubic Bezier segment described by the 2-tuples P0, ..., P3, return
    the y-value associated with the given x-value.

    Ex: getYfromXforBezSegment((10,0), (5,-5), (5,5), (0,0), 3.2)
    '''
    seconds_per_day = 24*60*60

    # Check if the second of the day is provided.
    # If provided, calculate y of the Bezier curve with that x
    # Otherwise, use the current second of the day
    if second_of_day is None:
        now = datetime.datetime.now()
        dt = datetime.timedelta(hours=now.hour, minutes=now.minute, seconds=now.second)
        seconds = dt.total_seconds()
    else:
        seconds = second_of_day

    # Shift the entire graph using 0 - 360 to determine the degree
    if shift_angle:
        percent_angle = shift_angle/360
        angle_seconds = percent_angle*seconds_per_day
        if seconds+angle_seconds > seconds_per_day:
            seconds_shifted = seconds+angle_seconds-seconds_per_day
        else:
            seconds_shifted = seconds+angle_seconds
        percent_of_day = seconds_shifted/seconds_per_day
    else:
        percent_of_day = seconds/seconds_per_day

    x = percent_of_day*(P0[0]-P3[0])

    # First, get the t-value associated with x-value, where t is the
    # parameterization of the Bezier curve and ranges from 0 to 1.
    # We need the coefficients of the polynomial describing cubic Bezier
    # (cubic polynomial in t)
    coefficients = [-P0[0] + 3*P1[0] - 3*P2[0] + P3[0],
                    3*P0[0] - 6*P1[0] + 3*P2[0],
                    -3*P0[0] + 3*P1[0],
                    P0[0] - x]
    # Find roots of the polynomial to determine the parameter t
    roots = np.roots(coefficients)
    # Find the root which is between 0 and 1, and is also real
    correct_root = None
    for root in roots:
        if np.isreal(root) and 0 <= root <= 1:
            correct_root = root
    # Check a valid root was found
    if correct_root is None:
        print('Error, no valid root found. Are you sure your Bezier curve '
              'represents a valid function when projected into the xy-plane?')
    param_t = correct_root
    # From the value for the t parameter, find the corresponding y-value
    # using the formula for cubic Bezier curves
    y = (1-param_t)**3*P0[1] + 3*(1-param_t)**2*param_t*P1[1] + 3*(1-param_t)*param_t**2*P2[1] + param_t**3*P3[1]
    assert np.isreal(y)
    # Typecast y from np.complex128 to float64   
    y = y.real
    return y

def kurt(x, opt):
    if opt=='kurt2':
        if np.all(x==0):
            K=0
            E=0
            return K
        x = x - np.mean(x)
        E = np.mean(np.abs(x)**2)
        if E < eps:
            K=0
            return K
        K = np.mean(np.abs(x)**4)/E**2
        if np.all(np.isreal(x)):
            K = K - 3
        else:
            K = K - 2
    if opt=='kurt1':
        if np.all(x==0):
            K=0
            E=0
            return K
        x = x - np.mean(x)
        E = np.mean(np.abs(x))
        if E < eps:
            K=0
            return K
        K = np.mean(np.abs(x)**2)/E**2
        if np.all(np.isreal(x)):
            K=K-1.57
        else:
            K=K-1.27
    return K

def read_matrix_sparse(filename, dtype=float, comments='#'):
    coo=np.loadtxt(filename, comments=comments, dtype=dtype)
    
    """Check if coo is (M, 3) ndarray"""
    if len(coo.shape)==2 and coo.shape[1]==3:
        row=coo[:, 0]
        col=coo[:, 1]
        values=coo[:, 2]

        """Check if imaginary part of row and col is zero"""
        if np.all(np.isreal(row)) and np.all(np.isreal(col)):
            row=row.real
            col=col.real

            """Check if first and second column contain only integer entries"""
            if np.all(is_integer(row)) and np.all(is_integer(col)):           

                """Convert row and col to int"""
                row=row.astype(int)
                col=col.astype(int)

                """Create coo-matrix"""
                A=scipy.sparse.coo_matrix((values,(row, col)))
                return A
            else:
                raise ValueError('File contains non-integer entries for row and col.')        
        else:
            raise ValueError('File contains complex entries for row and col.')
    else:
        raise ValueError('Given file is not a sparse matrix in coo-format.')

    def idct(self, X):
        '''Compute inverse discrete cosine transform of a 1-d array X'''
        
        N = len(X)
        w = np.sqrt(2*N)*np.exp(1j*np.arange(N)*np.pi/(2.*N))

        if (N%2==1) or (any(np.isreal(X))== False):
            w[0] = w[0] * np.sqrt(2)
            yy = np.zeros(2*N)
            yy[0:N] = w * X
            yy[N+1:2*N] = -1j * w[1:N] * X[1:N][::-1]
            y = fftpack.ifft(yy)
            x = y[0:N]
        else:
            w[0] = w[0]/np.sqrt(2)
            yy = X *w
            y = fftpack.ifft(yy)
            x = np.zeros(N)
            x[0:N:2] = y[0:N/2]
            x[1:N:2] = y[N:(N/2)-1:-1]

        if all(np.isreal(x)):
            x = np.real(x)

        return x

    def dct(self, x):
        '''Compute discrete cosine transform of 1-d array x'''

        #probably don't need this here anymore since it is in fftpack now    

        N = len(x)
        #calculate DCT weights
        w = (np.exp(-1j*np.arange(N)*np.pi/(2*N))/np.sqrt(2*N))
        w[0] = w[0]/np.sqrt(2)

        #test for odd or even function 
        if (N%2==1) or (any(np.isreal(x)) == False):
            y = np.zeros(2*N)
            y[0:N] = x
            y[N:2*N] = x[::-1]
            yy = fftpack.fft(y)
            yy = yy[0:N]
        else:
            y = np.r_[(x[0:N:2], x[N:0:-2])]
            yy = fftpack.fft(y)
            w = 2*w

        #apply weights
        X = w * yy

        if all(np.isreal(x)):
            X = np.real(X)

        return X

 def update_data(self, attr, old, new):
     mat = matrix(self.a_slider.value, self.b_slider.value,
                  self.c_slider.value, self.d_slider.value)
     evals, evecs = np.linalg.eig(mat)
     ev0, ev1 = evals
     evec0 = evecs.T[0]
     evec1 = evecs.T[1]
     trans0 = mat @ evec0
     trans1 = mat @ evec1
     self.eigen0.text = "eigen0 = " + \
         str(ev0) if np.isreal(ev0) else "eigen0 = None"
     self.eigen1.text = "eigen1 = " + \
         str(ev1) if np.isreal(ev1) else "eigen1 = None"
     transXs, transYs = mat @ np.array([self.Xs, self.Ys])
     self.source.data = dict(Xs=self.Xs, Ys=self.Ys,
                             transXs=transXs, transYs=transYs)
     data0x = [0, evec0[0]] if np.isreal(ev0) else [0, 0]
     data0y = [0, evec0[1]] if np.isreal(ev0) else [0, 0]
     data1x = [0, evec1[0]] if np.isreal(ev1) else [0, 0]
     data1y = [0, evec1[1]] if np.isreal(ev1) else [0, 0]
     trans0x = [0, trans0[0]] if np.isreal(ev0) else [0, 0]
     trans0y = [0, trans0[1]] if np.isreal(ev0) else [0, 0]
     trans1x = [0, trans1[0]] if np.isreal(ev1) else [0, 0]
     trans1y = [0, trans1[1]] if np.isreal(ev1) else [0, 0]
     self.evectors.data = dict(ev0x=data0x, ev0y=data0y,
                               ev1x=data1x, ev1y=data1y,
                               transev0x=trans0x, transev0y=trans0y,
                               transev1x=trans1x, transev1y=trans1y)

def sp_bin_by_space_and_time(data, frameTime, xybin=8, counterTime=40.0e-9):
    """Takes a dataset collected by the SSI photon counting fast camera and bins
    the data in time and xy.

    arguments:
        data - data to bin.  This is an (N, 4) array where N is the number of
            elements.
        frameTime - Amout of time between bins.  Measured in seconds
        xybin - amount to bin in x and y.  defaults to 8, which is a 512x512
            output.
        counterTime - clock time of the camera. Anton Tremsin says its 40 ns.

    returns:
        binnedData - a 3-dimension array (frames, y, x) of the binned data.
    """
    assert isinstance(data, np.ndarray), "data must be an ndarray"
    assert data.shape[1] == 4, "Second dimension must be 4."
    assert np.isreal(frameTime), "frameTime (%s) must be real" % repr(frameTime)
    assert np.isreal(counterTime), "counterTime (%s) must be real" % repr(counterTime)
    # TODO: This could be accelerated if written for a GPU

    firsttime = data[:,3].min()
    lasttime = data[:, 3].max()
    datalen = len(data[:, 0])
    
    nbins = int(np.ceil(float(lasttime - firsttime)*counterTime/frameTime))
    xybinnedDim = int(np.ceil(SP_MAX_SIZE/xybin))

    binnedData = np.zeros((nbins, xybinnedDim, xybinnedDim), dtype='int')
    bin_edges = np.linspace(firsttime, lasttime, nbins+1)
    # slightly increase the last bin.  For some reason the last bin is slightly too small, and the photons on the edges are not included
    bin_edges[-1] = bin_edges[-1]*1.01

    # If the frameTime is larger than the total acq time, then sum up everything
    if frameTime >= (lasttime-firsttime)*counterTime:
        histfn = lambda a,b: (len(a), b)
    else:
        # this is the default; bin data using histogram
        histfn = np.histogram

    for i in range(xybinnedDim):
        x = data[:, 0]
        y = data[:, 1]
        xc = np.argwhere( (i*xybin <= x) & (x < (i+1)*xybin) )
        idx_to_delete = np.array([])
        for j in range(xybinnedDim):
            yc = np.argwhere( (j*xybin <= y) & (y < (j+1)*xybin) )
            isect, a, b = intersect(xc, yc, assume_unique=True)
            if len(isect) != 0:
                binned, bin_edges = histfn(data[isect,3]+0.5, bin_edges)
                binnedData[:, i, j] = binned
            idx_to_delete = np.append(idx_to_delete, isect)
        data = np.delete(data, idx_to_delete, axis=0)

    if binnedData.sum() != datalen:
        print "warning: # of binned data photons (%d) do not match original dataset (%d)" % (binnedData.sum(), datalen)

    return binnedData

 def test_it(self):
     self.assertTrue(np.isreal(ef.return_real_part(1)))
     self.assertTrue(np.isreal(ef.return_real_part(1 + 0j)))
     self.assertTrue(np.isreal(ef.return_real_part(1 + 1e-20j)))
     self.assertRaises(TypeError, ef.return_real_part, None)
     self.assertRaises(TypeError, ef.return_real_part, (1, 2.0, 2 + 2j))
     self.assertRaises(TypeError, ef.return_real_part, [None, 2.0, 2 + 2j])
     self.assertRaises(ValueError, ef.return_real_part, [1, 2.0, 2 + 2j])
     self.assertRaises(ValueError, ef.return_real_part, 1 + 1e-10j)
     self.assertRaises(ValueError, ef.return_real_part, 1j)

def symmNormalization(MPS, chir, chic):
    omega, R = getLargestW(MPS, 'R')
    R = fixPhase(R)
    if np.isreal(R).all(): omega, R = omega.real, R.real
    print "wR", omega, np.isreal(R).all(), "R\n", R

    assym = np.linalg.norm(R - R.T.conj())
    print "assym R", assym

    Rvals, Rvecs = spla.eigh(R)
    Rvals_s = np.sqrt(abs(Rvals))
    Rvals_si = 1. / Rvals_s
    print "Rvals", Rvals

    Rs = np.dot(Rvecs, np.dot(np.diag(Rvals_s), Rvecs.T.conj()))
    ARs = np.tensordot(MPS, Rs, axes=([1,0]))
    Rsi = np.dot(Rvecs, np.dot(np.diag(Rvals_si), Rvecs.T.conj()))
    A1 = np.tensordot(Rsi, ARs, axes=([1,0]))
    A1 = np.transpose(A1, (0, 2, 1))

    omega, L = getLargestW(A1, 'L')
    L = fixPhase(L)
    if np.isreal(L).all(): omega, L = omega.real, L.real
    print "wL", omega, np.isreal(L).all(), "L\n", L

    assym = np.linalg.norm(L - L.T.conj())
    print "assym L", assym

    Lambda2, U = spla.eigh(L)
    A1U = np.tensordot(A1, U, axes=([1,0]))
    A2 = np.tensordot(U.T.conj(), A1U, axes=([1,0]))
    A2 = np.transpose(A2, (0, 2, 1))
    print "Lambda**2", Lambda2#, "\n", U

    Lambda = np.sqrt(abs(Lambda2))
    Lambdas = np.sqrt(Lambda)
    Lambdasi = 1. / Lambdas
    A2Lsi = np.tensordot(A2, np.diag(Lambdasi), axes=([1,0]))
    A3 = np.tensordot(np.diag(Lambdas), A2Lsi, axes=([1,0]))
    A3 = np.transpose(A3, (0, 2, 1))
    print "Lambda", Lambda

    nMPS = A3 / np.sqrt(omega)
    RealLambda = fixPhase(np.diag(Lambda))
    #print "nMPS", nMPS.shape, "\n", nMPS
    print "RealLambda", RealLambda.shape, "\n", RealLambda

    ######### CHECKING RESULT #########
    Trace = np.linalg.norm(RealLambda)
    ELambda = linearOpForR(nMPS, RealLambda).reshape(chir, chic)
    LambdaE = linearOpForL(nMPS, RealLambda).reshape(chir, chic)
    print "Trace(RealLambda)", Trace, "\nE|r)\n", ELambda, "\n(l|E\n", LambdaE

    return RealLambda, nMPS

 def _finalize_limits(self, axis, view, subplots, ranges):
     # Extents
     extents = self.get_extents(view, ranges)
     if extents and not self.overlaid:
         coords = [coord if np.isreal(coord) else np.NaN for coord in extents]
         if isinstance(view, Element3D) or self.projection == '3d':
             l, b, zmin, r, t, zmax = coords
             if not np.NaN in (zmin, zmax) and not zmin==zmax: axis.set_zlim((zmin, zmax))
         else:
             l, b, r, t = [coord if np.isreal(coord) else np.NaN for coord in extents]
         if not np.NaN in (l, r) and not l==r: axis.set_xlim((l, r))
         if not np.NaN in (b, t) and not b==t: axis.set_ylim((b, t))

def radialUndistort(rpp, k, quot=False, der=False):
    '''dqD
    takes distorted radius and returns the radius undistorted
    optioally it returns the undistortion quotient rp = rpp * q
    '''
    # polynomial coeffs, grade 7
    # # (k1,k2,p1,p2[,k3[,k4,k5,k6[,s1,s2,s3,s4[,τx,τy]]]])
    # polynomial coeffs, grade 7
    # # (k1,k2,p1,p2[,k3[,k4,k5,k6[,s1,s2,s3,s4[,τx,τy]]]])
    k.shape = -1
    poly = [[k[4], # k3
             0,
             k[1], # k2
             0,
             k[0], # k1
             0,
             1,
             -r] for r in rpp]
    
    # calculate roots
    rootsPoly = array([roots(p) for p in poly])
    
    # return flag, True if there is a suitable (real AND positive) solution
    rPRB = isreal(rootsPoly) & (0 <= real(rootsPoly))  # real Positive Real Bool
    retVal = any(rPRB, axis=1)
    
    rp = empty_like(rpp)
    
    if any(~ retVal): # if at least one case of non solution
        # calculate extrema of polyniomial
        rExtrema = roots([7*k[4], 0, 5*k[1], 0, 3*k[0], 0, 1])
        # select real extrema in positive side, keep smallest
        rRealPos = min(rExtrema.real[isreal(rExtrema) & (0<=rExtrema.real)])
        # assign to problematic values
        rp[~retVal] = rRealPos
    
    # choose minimum positive roots
    rp[retVal] = [min(rootsPoly[i, rPRB[i]].real)
                   for i in arange(rpp.shape[0])[retVal]]
    
    if der:
        # derivada de la directa
        q, dQdP, dQdK = radialDistort(rp, k, der=True)

        if quot:
            return q, retVal, dQdP, dQdK
        else:
            return rp, retVal, dQdP, dQdK
    else:
        if quot:
            return rp / rpp, retVal
        else:
            return rp, retVal

def kurt(x, opt):
    """
    Calculates kurtosis of a signal according to the option chosen

    :param x: signal
    :param opt: ['kurt1' | 'kurt2']
        * 'kurt1' = variance of the envelope magnitude
        * 'kurt2' = kurtosis of the complex envelope

    :type x: numpy array
    :type opt: string

    :rtype: float
    :returns: Kurtosis
    """
    if opt == 'kurt2':
        if np.all(x == 0):
            K = 0
            E = 0
            return K
        x = x - np.mean(x)
        E = np.mean(np.abs(x)**2)
        if E < eps:
            K = 0
            return K

        K = np.mean(np.abs(x)**4)/E**2

        if np.all(np.isreal(x)):
            K = K - 3
        else:
            K = K - 2

    if opt == 'kurt1':
        if np.all(x == 0):
            K = 0
            E = 0
            return K
        x = x - np.mean(x)
        E = np.mean(np.abs(x))
        if E < eps:
            K = 0
            return K

        K = np.mean(np.abs(x)**2)/E**2

        if np.all(np.isreal(x)):
            K = K-1.57
        else:
            K = K-1.27

    return K

def radialUndistort(rpp, k, quot=False, der=False):
    '''
    takes distorted radius and returns the radius undistorted
    optionally it returns the distortion quotioent rpp = rp * q
    '''
    # polynomial coeffs, grade 7
    # # (k1,k2,p1,p2[,k3[,k4,k5,k6[,s1,s2,s3,s4[,τx,τy]]]])
    # polynomial coeffs, grade 7
    # # (k1,k2,p1,p2[,k3[,k4,k5,k6[,s1,s2,s3,s4[,τx,τy]]]])
    
    k.shape = -1
    poly = [[k[3], 0, k[2], 0, k[1], 0, k[0], 0, 1, -r] for r in rpp]
    
    # calculate roots
    rootsPoly = array([roots(p) for p in poly])
    
    # return flag, True if there is a suitable (real AND positive) solution
    rPRB = isreal(rootsPoly) & (0 <= real(rootsPoly))  # real Positive Real Bool
    retVal = any(rPRB, axis=1)
    
    thetap = empty_like(rpp)
    
    if any(~retVal): # if at least one case of non solution
        # calculate extrema of polyniomial
        thExtrema = roots([9*k[3], 0, 7*k[2], 0, 5*k[1], 0, 3*k[0], 0, 1])
        # select real extrema in positive side, keep smallest
        thRealPos = min(thExtrema.real[isreal(thExtrema) & (0<=thExtrema.real)])
        # assign to problematic values
        thetap[~retVal] = thRealPos
    
    # choose minimum positive roots
    thetap[retVal] = [min(rootsPoly[i, rPRB[i]].real)
                   for i in arange(rpp.shape[0])[retVal]]
    
    rp = abs(tan(thetap))  # correct negative values
    # if theta angle is greater than pi/2, retVal=False
    retVal[thetap >= pi/2] = False
    
    if der:
        # derivada de la directa
        q, dQdP, dQdK = radialDistort(rp, k, quot, der)
        
        if quot:
            return q, retVal, dQdP, dQdK
        else:
            return rp, retVal, dQdP, dQdK
    else:
        if quot:
            # returns q
            return rpp / rp, retVal
        else:
            return rp, retVal

 def computeStationaryCovariance(self, ls):
     '''Compute the asymptotic covariance for a given linear system.
     (Eq. 104-106)
     '''
     A, H, G, Q, M, R = map(np.mat, [ls.getA(), ls.getH(), ls.getG(),
                                     ls.getQ(), ls.getM(), ls.getR()])
     try:
         Pprd = np.mat(dare(A.T, H.T, G * Q * G.T, M * R * M.T))
     except:
         raise AssertionError("Dare Unsolvable: "
                       + "The given system state is not stable/observable.")
     assert np.all(np.isreal(Pprd)) and np.all(Pprd == Pprd.T)
     Pest = Pprd - (Pprd * H.T) * (H * Pprd * H.T + M * R * M.T).I * (Pprd * H.T).T
     assert np.all(np.isreal(Pest)) and np.all(Pest == Pest.T)
     return Pest

    def multiply(self,x,mode):
        op1 = self.children[0]
        if np.isreal(x).all():
            # Purely real
            y = np.real(op1.applyMultiply(x,mode))
        elif np.isreal(1j*x).all():
            # Purely imaginary
            y = np.real(op1.applyMultiply(np.imag(x),mode)) * 1j
        else:
            # Mixed
            y = np.real(op1.applyMultiply(np.real(x),mode)) + np.real(op1.applyMultiply(np.imag(x),mode)) * 1j
            
        return y