def test_dtype_agreement():
    dtypes = [np.complex64, np.complex128, np.float32, np.float64]
    for dtype1 in dtypes:
        for dtype2 in dtypes:
            print dtype1,dtype2
            y_data = np.random.randn(10)+1j*np.random.randn(10)
            errors = np.random.randn(10)+1j*np.random.randn(10)
            with warnings.catch_warnings():
                warnings.simplefilter('ignore', np.ComplexWarning)
                y_data = y_data.astype(dtype1)
                errors = errors.astype(dtype2)
            x_data = np.linspace(100,110,10)
            if np.iscomplexobj(y_data) and np.iscomplexobj(errors):
                kid_readout.analysis.fitter.Fitter(x_data=x_data,y_data=y_data,errors=errors,
                                                   model=complex_dummy_model,
                                                   guess=complex_dummy_guess)
            elif (not np.iscomplexobj(y_data)) and (not np.iscomplexobj(errors)):
                kid_readout.analysis.fitter.Fitter(x_data=x_data,y_data=y_data,errors=errors,
                                                   model=kid_readout.analysis.fitter.line_model,
                                                   guess=kid_readout.analysis.fitter.line_guess)
            else:
                try:
                    kid_readout.analysis.fitter.Fitter(x_data=x_data,y_data=y_data,errors=errors,
                                                   model=complex_dummy_model,
                                                   guess=complex_dummy_guess)
                except TypeError:
                    pass

def map_coordinates(input, coordinates, output_type = None, output = None,
                order = 3, mode = 'constant', cval = 0.0, prefilter = True):
    """Apply an arbritrary coordinate transformation.

    The array of coordinates is used to find, for each point in the output,
    the corresponding coordinates in the input. The value of the input at
    that coordinates is determined by spline interpolation of the
    requested order.

    The shape of the output is derived from that of the coordinate
    array by dropping the first axis. The values of the array along
    the first axis are the coordinates in the input array at which the
    output value is found.  For example, if the input has dimensions
    (100,200,3), then the shape of coordinates will be (3,100,200,3),
    where coordinates[:,1,2,3] specify the input coordinate at which
    output[1,2,3] is found.

    Points outside the boundaries of the input are filled according to
    the given mode ('constant', 'nearest', 'reflect' or 'wrap'). The
    parameter prefilter determines if the input is pre-filtered before
    interpolation (necessary for spline interpolation of order >
    1). If False it is assumed that the input is already filtered.

    Example usage:
      >>> a = arange(12.).reshape((4,3))
      >>> print a
      [[  0.   1.   2.]
       [  3.   4.   5.]
       [  6.   7.   8.]
       [  9.  10.  11.]]
      >>> output = map_coordinates(a,[[0.5, 2], [0.5, 1]],order=1)
      >>> print output
      [ 2. 7.]

      Here, the interpolated value of a[0.5,0.5] gives output[0], while
      a[2,1] is output[1].
    """
    if order < 0 or order > 5:
        raise RuntimeError, 'spline order not supported'
    input = numpy.asarray(input)
    if numpy.iscomplexobj(input):
        raise TypeError, 'Complex type not supported'
    coordinates = numpy.asarray(coordinates)
    if numpy.iscomplexobj(coordinates):
        raise TypeError, 'Complex type not supported'
    output_shape = coordinates.shape[1:]
    if input.ndim < 1 or len(output_shape) < 1:
        raise RuntimeError, 'input and output rank must be > 0'
    if coordinates.shape[0] != input.ndim:
        raise RuntimeError, 'invalid shape for coordinate array'
    mode = _extend_mode_to_code(mode)
    if prefilter and order > 1:
        filtered = spline_filter(input, order, output = numpy.float64)
    else:
        filtered = input
    output, return_value = _ni_support._get_output(output, input,
                                        output_type, shape = output_shape)
    _nd_image.geometric_transform(filtered, None, coordinates, None, None,
               output, order, mode, cval, None, None)
    return return_value

def assert_array_almost_equal_nulp(x, y, nulp=1):
    """Compare two arrays relatively to their spacing. It is a relatively
    robust method to compare two arrays whose amplitude is variable.

    Note
    ----
    An assertion is raised if the following condition is not met:

        abs(x - y) <= nulps * spacing(max(abs(x), abs(y)))

    Parameters
    ----------
    x: array_like
        first input array
    y: array_like
        second input array
    nulp: int
        max number of unit in the last place for tolerance (see Note)
    """
    import numpy as np
    ax = np.abs(x)
    ay = np.abs(y)
    ref = nulp * np.spacing(np.where(ax > ay, ax, ay))
    if not np.all(np.abs(x-y) <= ref):
        if np.iscomplexobj(x) or np.iscomplexobj(y):
            msg = "X and Y are not equal to %d ULP" % nulp
        else:
            max_nulp = np.max(nulp_diff(x, y))
            msg = "X and Y are not equal to %d ULP (max is %g)" % (nulp, max_nulp)
        raise AssertionError(msg)

def write_memory_to_file(A, filename, mode='w', title='test'):
    """
    write memory to a h5 file
    h5 file contains root.real and root.imag(if A complex)
    best for transfer data with Matlab

    A: a ndarray, GPUArray or PitchArray
    filename: name of file to store
    mode: 'w' to start a new file
          'a' to append, leading dimension of A must be the
            same as the existing file

    file can be read by read_file or in matlab using h5read.m
    """
    h5file = tables.openFile(filename, mode, title)

    if (A.dtype == np.float32) or (A.dtype == np.complex64):
        tb = tables.Float32Atom
    elif (A.dtype == np.float64) or (A.dtype == np.complex128):
        tb = tables.Float64Atom
    elif A.dtype == np.int32:
        tb = tables.Int32Atom
    elif A.dtype == np.int64:
        tb = tables.Int64Atom
    else:
        TypeError("Write file error: unkown input dtype")

    if PYCUDA:
        if A.__class__.__name__ in ["GPUArray", "PitchArray"]:
            B = A.get()
        elif A.__class__.__name__ == "ndarray":
            B = A
        else:
            raise TypeError("Write file error: unkown input")
    else:
        if A.__class__.__name__ == "ndarray":
            B = A
        else:
            raise TypeError("Write file error: unkown input")

    shape = list(B.shape)
    shape[0] = 0

    if mode == 'w':
        if np.iscomplexobj(B):
            h5file.createEArray("/", "real", tb(), tuple(shape))
            h5file.createEArray("/", "imag", tb(), tuple(shape))
        else:
            h5file.createEArray("/", "real", tb(), tuple(shape))

    if np.iscomplexobj(B):
        h5file.root.real.append(B.real)
        h5file.root.imag.append(B.imag)
    else:
        h5file.root.real.append(B)

    h5file.close()

    if mode == 'w':
        print "file %s created" % (filename)

    def filter(self, array, *args, **kwargs):

        # Processed bandwith in percentages
        #------------------------------------------------------------------------
        sys = kwargs["meta"]        
        
        bw_proc_az = 1.33 *  sys['v0']     / sys['res_az']
        bw_proc_rg = 1.33 * (sys['c0']/2.) / sys['res_rg']
        
        percentage_az = 100.* bw_proc_az / (sys['prf']/sys['pre_az'])
        percentage_rg = 100.* bw_proc_rg /  sys['rsf']
        
        bw  = [percentage_az, percentage_rg]
        print("  bw=["+str(bw[0])+","+str(bw[1])+"] ... ")
        #------------------------------------------------------------------------        

        if array.ndim == 2 and np.iscomplexobj(array):
            return self.unweight2d(array, self.ov, bw)
        
        if array.ndim == 3 and np.iscomplexobj(array):
            p = array.shape
            for k in range(0,p[0]):
                array_temp = self.unweight2d(array[k,:,:], self.ov, bw)
                if k == 0:
                    s = array_temp.shape
                    array_new = np.empty((p[0],s[0],s[1]), dtype='complex64')
                array_new[k,:,:] = array_temp
            return array_new
        
        else:
            print("  ERROR: Bad input.")
            return None

    def fitToData(self, data):
        '''
        param data: numpy array where [:,0] is x and [:,1] is y
        '''
        x = data[:, 0][:, np.newaxis]
        y = data[:, 1][:, np.newaxis]
        D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
        S = np.dot(D.T, D)
        C = np.zeros([6, 6])
        C[0, 2] = C[2, 0] = 2; C[1, 1] = -1
        E, V = eig(np.dot(inv(S), C))
        n = np.argmax(np.abs(E))
        self.parameters = V[:, n]

        axes = self.ellipse_axis_length()
        self.a = axes[0]
        self.b = axes[1]
        self.angle = self.ellipse_angle_of_rotation()

        if not self.a or not self.b or self.parameters == None or np.iscomplexobj(self.parameters) or \
           math.isnan(self.a) or math.isnan(self.b) or math.isnan(self.ellipse_center()[0]) or \
           np.iscomplex(self.ellipse_center()[0]) or np.iscomplex(self.a) or np.iscomplex(self.b) or \
           np.iscomplexobj(self.angle):
            self.a = 0
            self.b = 0
            self.parameters = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
            self.angle = 0
            self.error = True

def subsample(args):
    """
    Rebin / Congrid variant
    """
    arr, shape, mode = args  # unpack arguments

    if mode == 'phase' and np.iscomplexobj(arr):
        arr = np.angle(arr)
    if np.iscomplexobj(arr):
        arr = np.abs(arr)

    if arr.shape == shape:
        return arr

    oshap = arr.shape

    for d in range(arr.ndim):
        n1 = shape[d]
        n2 = oshap[d]
        if n1 < n2:
            s = list(arr.shape)
            s.insert(d + 1, n2 // n1)
            s[d] = n1
            if mode == 'phase':
                arr = np.angle(np.exp(1j * arr.reshape(s)).mean(d + 1))
            elif mode == 'lables':
                arr = np.take(arr.reshape(s), 0, d + 1)
            else:
                arr = arr.reshape(s).mean(d + 1)
        else:
            arr = arr.repeat(n1 // n2, axis=d)
    return arr

    def _matvec(self, x):

        from bempp.api.utils.data_types import combined_type

        if x.ndim == 1:
            x_new = _np.expand_dims(x, 1)
            return self.matvec(x_new).ravel()

        if not self._fill_complete():
            raise ValueError("Not all rows or columns contain operators.")

        row_dim = 0
        res = _np.zeros((self.shape[0], x.shape[1]), dtype=combined_type(self.dtype, x.dtype))

        for i in range(self._m):
            col_dim = 0
            local_res = res[row_dim:row_dim + self._rows[i], :]
            for j in range(self._n):
                local_x = x[col_dim:col_dim + self._cols[j], :]
                if self._operators[i, j] is not None:
                    op_is_complex = _np.iscomplexobj(self._operators[i, j].dtype.type(1))
                    if _np.iscomplexobj(x) and not op_is_complex:
                        local_res[:] += (self._operators[i, j] * _np.real(local_x) +
                                         1j * self._operators[i, j] * _np.imag(local_x))
                    else:
                        local_res[:] += self._operators[i, j].dot(local_x)
                col_dim += self._cols[j]
            row_dim += self._rows[i]
        return res

 def __init__(self, f, data, model=default_model, guess=default_guess, functions=default_functions, 
              mask=None, errors=None, weight_by_errors=True):
     """
     Instantiate a resonator using our current best model.
     Parameter model is a function S_21(params, f) that returns the
     modeled values of S_21.
     Parameter guess is a function guess(f, data) that returns a
     good-enough initial guess at all of the fit parameters.
     Parameter functions is a dictionary that maps keys that are
     valid Python variables to functions that take a Parameters
     object as their only argument.
     Parameter mask is a boolean array of the same length as f and
     data; only points f[mask] and data[mask] are used to fit the
     data. The default is to use all data. Use this to exclude
     glitches or resonances other than the desired one.
     """
     if not np.iscomplexobj(data):
         raise TypeError("Resonator data should always be complex, but got real values")
     if errors is not None:
         if not np.iscomplexobj(errors):
             errors = errors*(1+1j)  # ensure errors is complex
     super(Resonator,self).__init__(f,data,model=model,guess=guess,functions=functions,mask=mask,
                                    errors=errors,weight_by_errors=weight_by_errors)
     if self.x_data.max() < 1e6:
         self.freq_units_MHz = True
     else:
         self.freq_units_MHz = False
     self.freq_data = self.x_data
     self.s21_data = self.y_data

def imshow(ax, x, y, z, *args, **kwargs):
    if (np.iscomplexobj(x)):
        x = np.asfarray(x.real)
    else:
        x = np.asfarray(x)
    if (np.iscomplexobj(y)):
        y = np.asfarray(y.real)
    else:
        y = np.asfarray(y)
    assert(len(x) == z.shape[1])
    assert(len(y) == z.shape[0])
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    if (np.iscomplexobj(z)):
        zabs = abs(z)
    else:
        zabs = z
    zabs = np.asfarray(zabs)
    # Use this to center pixel around (x,y) values
    extent = (x[0]-dx/2.0, x[-1]+dx/2.0, y[0]-dy/2.0, y[-1]+dy/2.0)
    # Use this to let (x,y) be the lower-left pixel location (upper-left when origin = 'lower' is not used)
    #extent = (x[0], x[-1], y[0], y[-1])
    im = ax.imshow(zabs, extent = extent, *args, **kwargs)
    imshow_show_z(ax, z, x, y)
    ax.set_xlim((x[0], x[-1]))
    ax.set_ylim((y[0], y[-1]))
    return im

    def __init__(self, freq, s21,
                 model=default_model, guess=default_guess, functions=default_functions,
                 mask=None, errors=None):
        """
        Fit a resonator using the given model.

        f: the frequencies used in a sweep.

        s21: the complex S_21 data taken at the given frequencies.

        model: a function S_21(params, f) that returns the modeled values of S_21.

        guess: a function guess(f, s21) that returns a good-enough initial guess at all of the fit parameters.

        functions: a dictionary that maps keys that are valid Python variables to functions that take a Parameters
        object as their only argument.

        mask: a boolean array of the same length as f and s21; only points f[mask] and s21[mask] are used to fit the
        data and the default is to use all data; use this to exclude glitches or resonances other than the desired one.
        """
        if not np.iscomplexobj(s21):
            raise TypeError("Resonator s21 must be complex.")
        if errors is not None and not np.iscomplexobj(errors):
            raise TypeError("Resonator s21 errors must be complex.")
        super(Resonator, self).__init__(freq, s21,
                                        model=model, guess=guess, functions=functions, mask=mask, errors=errors)
        self.freq_data = self.x_data
        self.s21_data = self.y_data
        self.freq_units_MHz = self.freq_data.max() < 1e6

    def get_jk_incore(self, cell=None, dm=None, hermi=1, verbose=logger.DEBUG, kpt=None):
        '''Get Coulomb (J) and exchange (K) following :func:`scf.hf.RHF.get_jk_`.

        *Incore* version of Coulomb and exchange build only.
        Currently RHF always uses PBC AO integrals (unlike RKS), since
        exchange is currently computed by building PBC AO integrals.
        '''
        if cell is None: cell = self.cell
        if dm is None: dm = self.make_rdm1()
        if kpt is None: kpt = self.kpt

        log = logger.Logger
        if isinstance(verbose, logger.Logger):
            log = verbose
        else:
            log = logger.Logger(cell.stdout, verbose)

        log.debug('JK PBC build: incore only with PBC integrals')

        if self._eri is None:
            log.debug('Building PBC AO integrals')
            if kpt is not None and pyscf.lib.norm(kpt) > 1.e-15:
                raise RuntimeError("Non-zero k points not implemented for exchange")
            self._eri = ao2mo.get_ao_eri(cell)

        if np.iscomplexobj(dm) or np.iscomplexobj(self._eri):
            vj, vk = dot_eri_dm_complex(self._eri, dm, hermi)
        else:
            vj, vk = pyscf.scf.hf.dot_eri_dm(self._eri, dm, hermi)

        return vj, vk

    def dot(self, coords_a, coords_b, frac_coords=False):
        """
        Compute the scalar product of vector(s).

        Args:
            coords_a, coords_b: Array-like objects with the coordinates.
            frac_coords (bool): Boolean stating whether the vector
                corresponds to fractional or cartesian coordinates.

        Returns:
            one-dimensional `numpy` array.
        """
        coords_a, coords_b = np.reshape(coords_a, (-1,3)), \
                             np.reshape(coords_b, (-1,3))

        if len(coords_a) != len(coords_b):
            raise ValueError("")

        if np.iscomplexobj(coords_a) or np.iscomplexobj(coords_b):
            raise TypeError("Complex array!")

        if not frac_coords:
            cart_a, cart_b = coords_a, coords_b
        else:
            cart_a = np.reshape([self.get_cartesian_coords(vec)
                                 for vec in coords_a], (-1,3))
            cart_b = np.reshape([self.get_cartesian_coords(vec)
                                 for vec in coords_b], (-1,3))

        return np.array([np.dot(a,b) for a,b in zip(cart_a, cart_b)])

def RCCSD(mf, frozen=0, mo_coeff=None, mo_occ=None):
    __doc__ = ccsd.CCSD.__doc__
    import numpy
    from pyscf import lib
    from pyscf import scf
    from pyscf.soscf import newton_ah
    from pyscf.cc import dfccsd

    if isinstance(mf, scf.uhf.UHF):
        raise RuntimeError('RCCSD cannot be used with UHF method.')
    elif isinstance(mf, scf.rohf.ROHF):
        lib.logger.warn(mf, 'RCCSD method does not support ROHF method. ROHF object '
                        'is converted to UHF object and UCCSD method is called.')
        return UCCSD(mf, frozen, mo_coeff, mo_occ)

    if isinstance(mf, newton_ah._CIAH_SOSCF) or not isinstance(mf, scf.hf.RHF):
        mf = scf.addons.convert_to_rhf(mf)

    if getattr(mf, 'with_df', None):
        return dfccsd.RCCSD(mf, frozen, mo_coeff, mo_occ)

    elif numpy.iscomplexobj(mo_coeff) or numpy.iscomplexobj(mf.mo_coeff):
        return rccsd.RCCSD(mf, frozen, mo_coeff, mo_occ)

    else:
        return ccsd.CCSD(mf, frozen, mo_coeff, mo_occ)

def test_c2c(comm):
    # this test requires pfft-python 0.1.16.

    pm = ParticleMesh(BoxSize=8.0, Nmesh=[8, 8], comm=comm, dtype='complex128')
    numpy.random.seed(1234)
    if comm.rank == 0:
        Npar = 100
    else:
        Npar = 0

    pos = 1.0 * (numpy.arange(Npar * len(pm.Nmesh))).reshape(-1, len(pm.Nmesh)) * (7, 7)
    pos %= (pm.Nmesh + 1)
    layout = pm.decompose(pos)

    npos = layout.exchange(pos)
    real = pm.paint(npos)

    complex = real.r2c()

    real2 = complex.c2r()

    assert numpy.iscomplexobj(real)
    assert numpy.iscomplexobj(real2)
    assert numpy.iscomplexobj(complex)

    assert_array_equal(complex.cshape, pm.Nmesh)
    assert_array_equal(real2.cshape, pm.Nmesh)
    assert_array_equal(real.cshape, pm.Nmesh)

    real.readout(npos)
    assert_almost_equal(numpy.asarray(real), numpy.asarray(real2), decimal=7)

    def __init__(self, x_data, y_data,
                 model=line_model, guess=line_guess, functions=default_functions,
                 mask=None, errors=None, method='leastsq', **minimize_keywords):
        """
        Arguments:

        model: a function y(params, x) that returns the modeled values.

        guess: a function guess(x_data, y_data) that returns a Parameters object containing an initial guess at the
        fit parameters.

        functions: a dictionary that maps keys that are valid Python variables to functions that take a Parameters
        object as their only argument.

        mask: a boolean array of the same length as f and data; only points x_data[mask] and y_data[mask] are used to
        fit the data, and the default is to use all data.

        errors: an array of the same size and data type as y_data with the corresponding error values;

        method: a string representing the fitting method for lmfit.minimize to use.

        minimize_keywords: keyword arguments that are passed directly to lmfit.minimize.

        Returns:

        A new Fitter using the given data and model.
        """
        self.x_data = x_data
        if np.iscomplexobj(y_data):
            y_data = y_data.astype('complex')  # promote data to complex128 if needed
            if errors is not None:
                if not np.iscomplexobj(errors):
                    raise TypeError(
                        "y_data and errors must both be complex or real, but got complex data with real errors.")
                errors = errors.astype('complex')
        else:  # data is real
            y_data = y_data.astype('float')  # promote data to float64 if needed
            if errors is not None:
                if np.iscomplexobj(errors):
                    raise TypeError(
                        "y_data and errors must both be complex or real, but got real data with complex errors.")
                errors = errors.astype('float')
        self.y_data = y_data
        self._model = model
        self._functions = functions
        self.method = method
        self.minimize_keywords = minimize_keywords
        if mask is None:
            self.mask = np.ones(x_data.shape, dtype=np.bool)
        else:
            self.mask = mask
        self.errors = errors
        if errors is None:
            self.residual = self._residual_without_errors
        else:
            self.residual = self._residual_with_errors
        self.fit(guess(x_data[self.mask], y_data[self.mask]))

def olafilt(b, x, zi=None):
    """
    Filter a one-dimensional array with an FIR filter

    Filter a data sequence, `x`, using a FIR filter given in `b`.
    Filtering uses the overlap-add method converting both `x` and `b`
    into frequency domain first.  The FFT size is determined as the
    next higher power of 2 of twice the length of `b`.

    Parameters
    ----------
    b : one-dimensional numpy array
        The impulse response of the filter
    x : one-dimensional numpy array
        Signal to be filtered
    zi : one-dimensional numpy array, optional
        Initial condition of the filter, but in reality just the
        runout of the previous computation.  If `zi` is None or not
        given, then zero initial state is assumed.

    Returns
    -------
    y : array
        The output of the digital filter.
    zf : array, optional
        If `zi` is None, this is not returned, otherwise, `zf` holds the
        final filter delay values.
    """

    L_I = b.shape[0]
    # Find power of 2 larger that 2*L_I (from abarnert on Stackoverflow)
    L_F = 2<<(L_I-1).bit_length()
    L_S = L_F - L_I + 1
    L_sig = x.shape[0]
    offsets = range(0, L_sig, L_S)

    # handle complex or real input
    if np.iscomplexobj(b) or np.iscomplexobj(x):
        fft_func = np.fft.fft
        ifft_func = np.fft.ifft
        res = np.zeros(L_sig+L_F, dtype=np.complex128)
    else:
        fft_func = np.fft.rfft
        ifft_func = np.fft.irfft
        res = np.zeros(L_sig+L_F)

    FDir = fft_func(b, n=L_F)

    # overlap and add
    for n in offsets:
        res[n:n+L_F] += ifft_func(fft_func(x[n:n+L_S], n=L_F)*FDir)

    if zi is not None:
        res[:zi.shape[0]] = res[:zi.shape[0]] + zi
        return res[:L_sig], res[L_sig:]
    else:
        return res[:L_sig]

    def compute(self):

        # make a copy for changes
        data = self.getData('in').copy()

        # convert complex to mag
        if np.iscomplexobj(data):
            data = np.abs(data)

        # get rid of negative numbers
        if data.min() < 0.:
            data -= data.min()

        # normalize the data
        data_min = data.min()
        data_max = data.max()
        data_range = data_max - data_min

        val = self.getAttr('L W F C:', 'val')

        new_min = data_range * val['floor'] / 100.0 + data_min
        new_max = data_range * val['ceiling'] / 100.0 + data_min

        data[data < new_min] = new_min
        data[data > new_max] = new_max

        data = data - math.fabs(new_min)
        data = data / (new_max - math.fabs(new_min)) * 255

        image = numpy2qimage(data.astype(np.uint8))
        if image.isNull():
            self.log.warn("Image Viewer: cannot load image")
        else:
            self.setAttr('Viewport:', val=image)

        line = self.getAttr('Viewport:', 'line')
        if line:
            a = self.getData('in')
            if np.iscomplexobj(a):
                a = np.abs(a)

            # Make a line with "l" points
            x0, y0 = line[0]
            x1, y1 = line[1]
            l = int(np.hypot(x1 - x0, y1 - y0))
            x, y = np.linspace(x0, x1, l), np.linspace(y0, y1, l)

            # Extract the values along the line, using cubic interpolation
            arr = ndimage.map_coordinates(a, np.vstack((x, y)), order=3)

            self.setAttr('Cross Section', val=[arr])
            self.setAttr('Cross Section', visible=True)
        else:
            self.setAttr('Cross Section', visible=False)

        return(0)

    def plot(self, *args, **kwargs):
        if "projection" in kwargs:
            projection = kwargs.pop("projection")
        else:
            projection = self.name
        vars = args[:2]
        args = args[2:]

        if len(vars) == 2 and isinstance(vars[1], (str, unicode)):
            args = (vars[1],) + args
            vars = vars[:1]

        if ((len(vars) == 1 and
             isinstance(vars[0], hfarray) and
             len(vars[0].dims) >= 1)):
            y = vars[0]
            x = hfarray(y.dims[0])
            vars = (x, y)
            if self.HFTOOLS_default_x_name is None:
                self.HFTOOLS_default_x_name = y.dims[0].name
                fmt = self.axes.xaxis.get_major_formatter()
                if hasattr(fmt, "update_template"):
                    fmt.default_label = self.HFTOOLS_default_x_name
                    fmt.update_template()

        if len(vars) == 1:
            y = vars[0]
            if projection in _projfun:
                x, y = _projfun[projection](None, y)
                return Axes.plot(self, y, *args, **kwargs)
            elif np.iscomplexobj(y):
                return Axes.plot(self, y.real, y.imag, *args, **kwargs)
            else:
                return Axes.plot(self, y, *args, **kwargs)
        elif len(vars) == 2:
            x = vars[0]
            y = vars[1]
            xunit = getattr(x, "unit", None)
            yunit = getattr(y, "unit", None)

            if projection in _projfun:
                x, y = _projfun[projection](x, y)
                lines = self._plot_helper(x, y, *args, **kwargs)
            elif np.iscomplexobj(y):
                xunit = yunit
                lines = self._plot_helper(y.real, y.imag, *args, **kwargs)
            else:
                lines = self._plot_helper(x, y, *args, **kwargs)
            if xunit:
                self.set_xlabel_unit(xunit)
            if yunit:
                self.set_ylabel_unit(yunit)
            return lines
        else:
            raise Exception("Missing plot data")

    def __test(y, truth, max_size, axis):

        ac = librosa.autocorrelate(y, max_size=max_size, axis=axis)

        my_slice = [slice(None)] * truth.ndim
        if max_size is not None and max_size <= y.shape[axis]:
            my_slice[axis] = slice(min(max_size, y.shape[axis]))

        if not np.iscomplexobj(y):
            assert not np.iscomplexobj(ac)

        assert np.allclose(ac, truth[my_slice])