def call(self, args, axis=0, out=None, chunksize=1024 * 1024, **kwargs):
        """ axis is the axis to chop it off.
            if self.altreduce is set, the results will
            be reduced with altreduce and returned
            otherwise will be saved to out, then return out.
        """
        if self.altreduce is not None:
            ret = [None]
        else:
            if out is None :
                if self.outdtype is not None:
                    dtype = self.outdtype
                else:
                    try:
                        dtype = numpy.result_type(*[args[i] for i in self.ins] * 2)
                    except:
                        dtype = None
                out = sharedmem.empty(
                        numpy.broadcast(*[args[i] for i in self.ins] * 2).shape,
                        dtype=dtype)
        if axis != 0:
            for i in self.ins:
                args[i] = numpy.rollaxis(args[i], axis)
            out = numpy.rollaxis(out, axis)
        size = numpy.max([len(args[i]) for i in self.ins])
        with sharedmem.MapReduce() as pool:
            def work(i):
                sl = slice(i, i+chunksize)
                myargs = args[:]
                for j in self.ins:
                    try: 
                        tmp = myargs[j][sl]
                        a, b, c = sl.indices(len(args[j]))
                        myargs[j] = tmp
                    except Exception as e:
                        print tmp
                        print j, e
                        pass
                if b == a: return None
                rt = self.ufunc(*myargs, **kwargs)
                if self.altreduce is not None:
                    return rt
                else:
                    out[sl] = rt
            def reduce(rt):
                if self.altreduce is None:
                    return
                if ret[0] is None:
                    ret[0] = rt
                elif rt is not None:
                    ret[0] = self.altreduce(ret[0], rt)

            pool.map(work, range(0, size, chunksize), reduce=reduce)

        if self.altreduce is None:
            if axis != 0:
                out = numpy.rollaxis(out, 0, axis + 1)
            return out                
        else:
            return ret[0]

def shifted_corr(reference, image, displacement):
    """Calculate the correlation between the reference and the image shifted
    by the given displacement.

    Parameters
    ----------
    reference : np.ndarray
    image : np.ndarray
    displacement : np.ndarray

    Returns
    -------
    correlation : float

    """

    ref_cuts = np.maximum(0, displacement)
    ref = reference[ref_cuts[0]:, ref_cuts[1]:, ref_cuts[2]:]
    im_cuts = np.maximum(0, -displacement)
    im = image[im_cuts[0]:, im_cuts[1]:, im_cuts[2]:]
    s = np.minimum(im.shape, ref.shape)
    ref = ref[:s[0], :s[1], :s[2]]
    im = im[:s[0], :s[1], :s[2]]
    ref -= nanmean(ref.reshape(-1, ref.shape[-1]), axis=0)
    ref = np.nan_to_num(ref)
    im -= nanmean(im.reshape(-1, im.shape[-1]), axis=0)
    im = np.nan_to_num(im)
    assert np.all(np.isfinite(ref)) and np.all(np.isfinite(im))
    corr = nanmean(
        [old_div(np.sum(i * r), np.sqrt(np.sum(i * i) * np.sum(r * r))) for
         i, r in zip(np.rollaxis(im, -1), np.rollaxis(ref, -1))])
    return corr

def write_raw(filename, signal, record_by):
    """Writes the raw file object

    Parameters:
    -----------
    filename : string
        the filename, either with the extension or without it
    record_by : string
     'vector' or 'image'

        """
    filename = os.path.splitext(filename)[0] + '.raw'
    dshape = signal.data.shape
    data = signal.data
    if len(dshape) == 3:
        if record_by == 'vector':
            np.rollaxis(
                data, signal.axes_manager._slicing_axes[0].index_in_array, 3
                        ).ravel().tofile(filename)
        elif record_by == 'image':
            data = np.rollaxis(
                data, signal.axes_manager._non_slicing_axes[0].index_in_array, 0
                        ).ravel().tofile(filename)
    elif len(dshape) == 2:
        if record_by == 'vector':
            np.rollaxis(
                data, signal.axes_manager._slicing_axes[0].index_in_array, 2
                        ).ravel().tofile(filename)
        elif record_by in ('image', 'dont-care'):
            data.ravel().tofile(filename)
    elif len(dshape) == 1:
        data.ravel().tofile(filename)

    def extract(self, X, batch_size = 500):
        assert self._z is not None, "Must be trained before calling extract"
        X_size, X_w, X_h, X_channel = X.shape
        total_dimension = np.prod(self._part_shape)
        
        coded_result = np.zeros((X_size, X_w - self._part_shape[0] + 1, X_h - self._part_shape[1] + 1, self._num_features))
        
        X_input = tf.placeholder("float32", [None, total_dimension + 1], name = "X_input")
        fc_1 = tf.matmul(X_input, tf.transpose(self._z, [1, 0], name = "transpose"))
        
        code_function = tf.exp(fc_1)
        
        for i in range(X_size // batch_size):
            end = min((i + 1) * batch_size, X_size)
            start = i * batch_size
            X_select = X[start: end]
            code_x = np.ones((end-start, coded_result.shape[1], coded_result.shape[2], total_dimension + 1))
            norm_x = np.zeros((end-start, coded_result.shape[1], coded_result.shape[2]))
            for m in range(coded_result.shape[1]):
                for n in range(coded_result.shape[2]):
                    selected_patches = X_select[:,m:m+self._part_shape[0], n:n+self._part_shape[1], :].reshape(-1, total_dimension)
                    patches_norm = np.sqrt(np.sum(selected_patches ** 2, axis = 1))
                    patches_norm = np.clip(patches_norm, a_min = 0.00001, a_max = 10)
                    code_x[:, m, n, :total_dimension] = np.array([selected_patches[k] / patches_norm[k] for k in range(end-start)])
                    norm_x[:, m, n] = patches_norm
            code_x = code_x.reshape(-1, total_dimension + 1)
            feed_dict = {}
            feed_dict[X_input] = code_x
            tmp_coded_result = self._sess.run(code_function, feed_dict = feed_dict)
            reshape_result = tmp_coded_result.reshape(end-start, coded_result.shape[1], coded_result.shape[2], self._num_features)
            coded_result[start:end] = np.rollaxis(np.rollaxis(reshape_result, 3, 0) / norm_x, 0, 4)
            print norm_x

        return coded_result

 def _run_interface(self, runtime):
     img = nb.load(self.inputs.in_file[0])
     header = img.get_header().copy()
     vollist = [nb.load(filename) for filename in self.inputs.in_file]
     data = np.concatenate([vol.get_data().reshape(
         vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3)
     if data.dtype.kind == 'i':
         header.set_data_dtype(np.float32)
         data = data.astype(np.float32)
     if isdefined(self.inputs.regress_poly):
         timepoints = img.get_shape()[-1]
         X = np.ones((timepoints, 1))
         for i in range(self.inputs.regress_poly):
             X = np.hstack((X, legendre(
                 i + 1)(np.linspace(-1, 1, timepoints))[:, None]))
         betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2))
         datahat = np.rollaxis(np.dot(X[:, 1:],
                                      np.rollaxis(
                                          betas[1:, :, :, :], 0, 3)),
                               0, 4)
         data = data - datahat
         img = nb.Nifti1Image(data, img.get_affine(), header)
         nb.save(img, self._gen_output_file_name('detrended'))
     meanimg = np.mean(data, axis=3)
     stddevimg = np.std(data, axis=3)
     tsnr = meanimg / stddevimg
     img = nb.Nifti1Image(tsnr, img.get_affine(), header)
     nb.save(img, self._gen_output_file_name())
     img = nb.Nifti1Image(meanimg, img.get_affine(), header)
     nb.save(img, self._gen_output_file_name('mean'))
     img = nb.Nifti1Image(stddevimg, img.get_affine(), header)
     nb.save(img, self._gen_output_file_name('stddev'))
     return runtime

def load_data(trainingData, trainingLabel,
              testingData, testingLabel,
              resize = False, size = 100, dataset = "IKEA_PAIR"):
    trainingData = os.environ[dataset] + trainingData
    trainingLabel = os.environ[dataset] + trainingLabel
    testingData = os.environ[dataset] + testingData
    testingLabel = os.environ[dataset] + testingLabel

    X_train = np.array(np.load(trainingData),
                       dtype = np.float32)
    Y_train = np.array(np.load(trainingLabel),
                       dtype = np.uint8)
    X_test = np.array(np.load(testingData),
                      dtype = np.float32)
    Y_test = np.array(np.load(testingLabel),
                      dtype = np.uint8)

    print("resizing....")
    if resize:
        X_train = np.array([misc.imresize(X_train[i],
                                          size = (size, size, 3)) /255.0
                            for i in range(X_train.shape[0])], dtype=np.float32)
        X_test = np.array([misc.imresize(X_test[i],
                                         size = (size, size, 3)) /255.0
                           for i in range(X_test.shape[0])], dtype=np.float32)
        np.save(trainingData + "_100.npy", X_train)
        np.save(testingData + "_100.npy", X_test)

    X_train = np.rollaxis(X_train, 3, 1)
    X_test = np.rollaxis(X_test, 3, 1)

    print("downresizing....")


    return X_train, Y_train, X_test, Y_test

def loadData(folder):
	data = dp.DataProcessor()
	#plotData(data)
	nsp = data.normalizedSequencesPerCell()
	'''Select last 5 genes as target values'''
	oSamples = np.array(nsp[:,1:,3:7], dtype='float32')
	'''Bring array into shape (n_steps, n_samples, n_genes)'''
	oSamples = np.rollaxis(oSamples, 0, 2)
	'''Select first time point of last 5 genes as initial condition'''
	c0Samples = np.array(nsp[:,0,3:7], dtype='float32')
	
	stepsPerUnit = 1/dt
	numUnits = oSamples.shape[1]
	totalSteps = numUnits*stepsPerUnit
	
	'''Create input genes array'''
	interpGenes = []
	for g in xrange(3):
		interpGenes.append(interpolateSingleGene(nsp[:,:,g], totalSteps))
	
	iSamples = np.array(interpGenes, dtype='float32')
	
	'''Bring array into shape (n_steps, n_samples, n_genes)'''
	iSamples = np.rollaxis(iSamples, 0, 3).clip(min = 0)
	iSamples = np.rollaxis(iSamples, 0, 2)
	
	np.save(folder+"data/simpleInputSequence.npy", iSamples)
	np.save(folder+"data/simpleOutputSequence.npy", oSamples)
	np.save(folder+"data/simpleStartSequence.npy", c0Samples)
	
	return iSamples, oSamples, c0Samples

def _to_rgb_uint8(image, autoscale):
    if autoscale is None:
        autoscale = image.dtype != np.uint8

    if autoscale:
        image = (normalize(image) * 255).astype(np.uint8)
    elif image.dtype != np.uint8:
        if np.issubdtype(image.dtype, np.integer):
            max_value = np.iinfo(image.dtype).max
            # sometimes 12-bit images are stored as unsigned 16-bit
            if max_value == 2**16 - 1 and image.max() < 2**12:
                max_value = 2**12 - 1
            image = (image / max_value * 255).astype(np.uint8)
        else:
            image = (image * 255).astype(np.uint8)

    ndim = image.ndim
    shape = image.shape
    if ndim == 3 and shape.count(3) == 1:
        # This is a color image. Ensure that the color axis is axis 2.
        color_axis = shape.index(3)
        image = np.rollaxis(image, color_axis, 3)
    elif image.ndim == 3 and shape.count(4) == 1:
        # This is an RGBA image. Ensure that the color axis is axis 2, and 
        # drop the A values.
        color_axis = shape.index(4)
        image = np.rollaxis(image, color_axis, 3)[:, :, :3]
    elif ndim == 2:
        # Expand into color to satisfy moviepy's expectation
        image = np.repeat(image[:, :, np.newaxis], 3, axis=2)
    else:
        raise ValueError("Images have the wrong shape.")

    return np.asarray(image)

def read_data_sets(data_dir, distortion=True, dtype=np.float32, training_num=18000):
    global NUM_EXAMPLES_PER_EPOCH_FOR_TRAIN
    NUM_EXAMPLES_PER_EPOCH_FOR_TRAIN = training_num

    train_image = np.array(np.load(os.path.join(data_dir, "real_background_10_class_train.npy")).reshape(-1, 32, 32, 3), dtype=dtype)
    train_image = np.rollaxis(train_image, 3, 1)
    train_labels = np.array(np.load(os.path.join(data_dir, "10_class_train_label.npy")),dtype=dtype)
    total_num_train = train_image.shape[0]
    random_index = np.arange(total_num_train)
    np.random.shuffle(random_index)
    train_image = train_image[random_index]
    train_labels = train_labels[random_index]
    
    train_image = train_image[:NUM_EXAMPLES_PER_EPOCH_FOR_TRAIN]
    train_labels = train_labels[:NUM_EXAMPLES_PER_EPOCH_FOR_TRAIN]
    print(train_image.shape, train_labels.shape)

    test_image = np.array(np.load(os.path.join(data_dir, "real_background_10_class_test.npy")).reshape(-1, 32, 32, 3), dtype=dtype)
    test_image = np.rollaxis(test_image, 3, 1)
    test_labels = np.array(np.load(os.path.join(data_dir, "10_class_test_label.npy")), dtype=dtype)

    print(test_image.shape, test_labels.shape)

    train = DataSet(train_image, train_labels, distortion=distortion)
    test = DataSet(test_image, test_labels, test=True)

    Datasets = collections.namedtuple('Datasets', ['train', 'test'])

    return Datasets(train = train, test = test)

 def vec_func(p, x):
     # Move last dimension to front (becomes sequence of column arrays)
     column_seq_x = np.rollaxis(np.asarray(x), -1)
     # Get corresponding sequence of output column arrays
     column_seq_y = np.array([func(p, xx) for xx in column_seq_x])
     # Move column dimension back to the end
     return np.rollaxis(column_seq_y, 0, len(column_seq_y.shape))

def time_subset_awot_dict(time, data, start_time, end_time, time_axis=0):
    '''
    Get the variable from the fields dictionary.
    Subset the time when in time series format.

    Parameters
     ----------
    time : dict
        AWOT time dictionary
    data : dict
        AWOT data dictionary.
    start_time : str
        UTC time to use as start time for subsetting in datetime format.
        (e.g. 2014-08-20 12:30:00)
    end_time : str
        UTC time to use as an end time for subsetting in datetime format.
        (e.g. 2014-08-20 16:30:00)
    '''
    # Check to see if time is subsetted
    dt_start = _get_start_datetime(time, start_time)
    dt_end = _get_end_datetime(time, end_time)
    datasub = data.copy()
    if time_axis > 0:
        np.rollaxis(datasub['data'], time_axis)
    datasub['data'] = data['data'][(time['data'] >= dt_start) &
                                   (time['data'] <= dt_end), ...]
    return datasub

def bspread(X, spread='box', radius=1, first_axis=False):
    """
    Spread binary edges.

    Parameters
    ----------
    X : ndarray  (3D or 4D)
        Binary edges to spread. Shape should be ``(rows, cols, A)`` or ``(N, rows, cols, A)``, where `A` is the number of edge features.
    first_axis: bool
         If True, the images will be assumed to be ``(A, rows, cols)`` or ``(N, A, rows, cols)``. 
    spread : 'box', 'orthogonal', None 
        If set to `'box'` and `radius` is set to 1, then an edge will appear if any of the 8 neighboring pixels detected an edge. This is equivalent to inflating the edges area with 1 pixel. The size of the box is dictated by `radius`. 
        If `'orthogonal'`, then the features will be extended by `radius` perpendicular to the direction of the edge feature (i.e. along the gradient).
    radius : int
        Controls the extent of the inflation, see above.
    """
    single = X.ndim == 3
    if single:
        X = X.reshape((1,) + X.shape) 
    if not first_axis:
        X = np.rollaxis(X, 3, start=1)

    Xnew = array_bspread(X, spread, radius)

    if not first_axis:
        Xnew = np.rollaxis(Xnew, 1, start=4)

    if single:
        Xnew = Xnew.reshape(Xnew.shape[1:]) 

    return Xnew

def get_dH2(lab1, lab2):
    """squared hue difference term occurring in deltaE_cmc and deltaE_ciede94

    Despite its name, "dH" is not a simple difference of hue values.  We avoid
    working directly with the hue value, since differencing angles is
    troublesome.  The hue term is usually written as:
        c1 = sqrt(a1**2 + b1**2)
        c2 = sqrt(a2**2 + b2**2)
        term = (a1-a2)**2 + (b1-b2)**2 - (c1-c2)**2
        dH = sqrt(term)

    However, this has poor roundoff properties when a or b is dominant.
    Instead, ab is a vector with elements a and b.  The same dH term can be
    re-written as:
        |ab1-ab2|**2 - (|ab1| - |ab2|)**2
    and then simplified to:
        2*|ab1|*|ab2| - 2*dot(ab1, ab2)
    """
    lab1 = np.asarray(lab1)
    lab2 = np.asarray(lab2)
    a1, b1 = np.rollaxis(lab1, -1)[1:3]
    a2, b2 = np.rollaxis(lab2, -1)[1:3]

    # magnitude of (a, b) is the chroma
    C1 = np.hypot(a1, b1)
    C2 = np.hypot(a2, b2)

    term = (C1 * C2) - (a1 * a2 + b1 * b2)
    return 2 * term

def deltaE_cie76(lab1, lab2):
    """Euclidean distance between two points in Lab color space

    Parameters
    ----------
    lab1 : array_like
        reference color (Lab colorspace)
    lab2 : array_like
        comparison color (Lab colorspace)

    Returns
    -------
    dE : array_like
        distance between colors `lab1` and `lab2`

    References
    ----------
    .. [1] https://en.wikipedia.org/wiki/Color_difference
    .. [2] A. R. Robertson, "The CIE 1976 color-difference formulae,"
           Color Res. Appl. 2, 7-11 (1977).
    """
    lab1 = np.asarray(lab1)
    lab2 = np.asarray(lab2)
    L1, a1, b1 = np.rollaxis(lab1, -1)[:3]
    L2, a2, b2 = np.rollaxis(lab2, -1)[:3]
    return np.sqrt((L2 - L1) ** 2 + (a2 - a1) ** 2 + (b2 - b1) ** 2)

def cross(a, b, axisa=0, axisb=0, axisc=0):
    """Vector cross product along specified axes of ndarrays."""
    if not hasattr(_np, 'einsum'):
        import nein
        return nein.cross(a, b, axisa, axisb, axisc)
    a, b = _asarray(a, b)
    if (a.ndim != b.ndim and
            a.shape not in [(2,), (3,)] and
            b.shape not in [(2,), (3,)]):
        return _np.cross(a, b, axisa=axisa, axisb=axisb, axisc=axisc)
    axisa, axisb = _normalize_indices(a, b, axisa, axisb)
    n = a.shape[axisa]
    if n not in [2, 3]:
        raise NotImplementedError(
            "Only 2D and 3D cross products are implemented")
    if n != b.shape[axisb]:
        raise ValueError(_error(a, b, axisa, axisb))
    if n == 2:
        return _cross2d(a, b, axisa, axisb, axisc)
    strb = '%sj%s' % (_LS[:axisa], _LS[axisa:len(a.shape) - 1])
    strc = 'ik%s' % strb.replace('j', '')
    a = _ein(eijk, a, 'ijk', strb, strc)
    stra = strc
    strb = '%sk%s' % (_LS[:axisb], _LS[axisb:len(b.shape) - 1])
    series = ''.join(sorted(x for x in _LS if x in stra or x in strb))
    if axisc < 0:
        axisc += len(series) + 1
    strc = '%si%s' % (series[:axisc], series[axisc:])
    mask = _collapse_mask(_np.rollaxis(a, axisa),
                          _collapse_mask(_np.rollaxis(b, axisb)))
    return _ein(a, b, stra, strb, strc, mask=mask, mask_axes=(axisc,))

def load_svhn_src(roll=True, batchsize=600):
    """
    TODO : read again and rework it !
    """
    data = io.loadmat(os.path.join(data_dir,'train_32x32.mat'))
    X = data['X']
    y = data['y']
    X = np.rollaxis(X, 3)
    if roll:
        X = np.rollaxis(X, 3, 1)

    s1 = 50000
    s2 = 60000
    X_train, X_val, X_test = X[:s1], X[s1:s2], X[s2:]
    y_train, y_val, y_test = y[:s1], y[s1:s2], y[s2:]

    data = {
            'X_train': X_train,
            'y_train': y_train,
            'X_val': X_val,
            'y_val': y_val,
            'X_test': X_test,
            'y_test': y_test,
            'batchsize': batchsize,
            }
    return data

def mtimesm(a, b, axisa=0, axisb=0, axisc=0,
            transposea=False, transposeb=False, transposec=False, **kwargs):
    """Matrix/matrix multiplication along specified axes of ndarrays."""
    if not hasattr(_np, 'einsum'):
        import nein
        return nein.mtimesm(a, b, axisa, axisb, axisc,
                            transposea, transposeb, transposec, **kwargs)
    a, b = _asarray(a, b)
    axisa, axisb = _normalize_indices(a, b, axisa, axisb)
    if a.shape[axisa + 1] != b.shape[axisb]:
        raise ValueError(_error(a, b, axisa, axisb))
    transposea = kwargs.get('ta', transposea)
    transposeb = kwargs.get('tb', transposeb)
    transposec = kwargs.get('tc', transposec)
    stra = '%s%s%s' % (_LS[:axisa], 'ji' if transposea else 'ij',
                       _LS[axisa:len(a.shape) - 2])
    strb = '%s%s%s' % (_LS[:axisb], 'kj' if transposeb else 'jk',
                       _LS[axisb:len(b.shape) - 2])
    series = ''.join(sorted(x for x in _LS if x in stra or x in strb))
    if axisc < 0:
        axisc += len(series) + 2
    strc = '%s%s%s' % (series[:axisc], 'ki' if transposec else 'ik',
                       series[axisc:])
    mask = None
    for x, ax in (a, axisa), (b, axisb):
        if isinstance(x, _np.ma.MaskedArray):
            if x.ndim > 2:
                x = _np.rollaxis(_np.rollaxis(a, ax), ax + 1, 1)
                for row in x:
                    mask = _collapse_mask(row, mask)
            elif x.mask.any():
                mask = True
            elif mask is None:
                mask = False
    return _ein(a, b, stra, strb, strc, mask=mask, mask_axes=(axisc, axisc + 1))

def icrs_to_cirs(icrs_coo, cirs_frame):
    # first set up the astrometry context for ICRS<->CIRS
    astrom, eo = erfa.apci13(*get_jd12(cirs_frame.obstime, 'tdb'))

    if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
        # if no distance, just do the infinite-distance/no parallax calculation
        usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
        i_ra = usrepr.lon.to(u.radian).value
        i_dec = usrepr.lat.to(u.radian).value
        cirs_ra, cirs_dec = atciqz(i_ra, i_dec, astrom)

        newrep = UnitSphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
                                             lon=u.Quantity(cirs_ra, u.radian, copy=False),
                                             copy=False)
    else:
        # When there is a distance,  we first offset for parallax to get the
        # astrometric coordinate direction and *then* run the ERFA transform for
        # no parallax/PM. This ensures reversiblity and is more sensible for
        # inside solar system objects
        newxyz = icrs_coo.cartesian.xyz
        newxyz = np.rollaxis(newxyz, 0, newxyz.ndim) - astrom['eb'] * u.au
        # roll xyz back to the first axis
        newxyz = np.rollaxis(newxyz, -1, 0)
        newcart = CartesianRepresentation(newxyz)

        srepr = newcart.represent_as(SphericalRepresentation)
        i_ra = srepr.lon.to(u.radian).value
        i_dec = srepr.lat.to(u.radian).value
        cirs_ra, cirs_dec = atciqz(i_ra, i_dec, astrom)

        newrep = SphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
                                         lon=u.Quantity(cirs_ra, u.radian, copy=False),
                                         distance=srepr.distance, copy=False)

    return cirs_frame.realize_frame(newrep)

 def setup_method(self, method):
     # Create three signals with dimensions:
     # s1 : <BaseSignal, title: , dimensions: (4, 3, 2|2, 3)>
     # s2 : <BaseSignal, title: , dimensions: (2, 3|4, 3, 2)>
     # s12 : <BaseSignal, title: , dimensions: (2, 3|4, 3, 2)>
     # Where s12 data is transposed in respect to s2
     dc1 = np.random.random((2, 3, 4, 3, 2))
     dc2 = np.rollaxis(np.rollaxis(dc1, -1), -1)
     s1 = signals.BaseSignal(dc1.copy())
     s2 = signals.BaseSignal(dc2)
     s12 = signals.BaseSignal(dc1.copy())
     for i, axis in enumerate(s1.axes_manager._axes):
         if i < 3:
             axis.navigate = True
         else:
             axis.navigate = False
     for i, axis in enumerate(s2.axes_manager._axes):
         if i < 2:
             axis.navigate = True
         else:
             axis.navigate = False
     for i, axis in enumerate(s12.axes_manager._axes):
         if i < 3:
             axis.navigate = False
         else:
             axis.navigate = True
     self.s1 = s1
     self.s2 = s2
     self.s12 = s12

def cirs_to_icrs(cirs_coo, icrs_frame):
    srepr = cirs_coo.represent_as(UnitSphericalRepresentation)
    cirs_ra = srepr.lon.to(u.radian).value
    cirs_dec = srepr.lat.to(u.radian).value

    # set up the astrometry context for ICRS<->cirs and then convert to
    # astrometric coordinate direction
    astrom, eo = erfa.apci13(*get_jd12(cirs_coo.obstime, 'tdb'))
    i_ra, i_dec = aticq(cirs_ra, cirs_dec, astrom)

    if cirs_coo.data.get_name() == 'unitspherical' or cirs_coo.data.to_cartesian().x.unit == u.one:
        # if no distance, just use the coordinate direction to yield the
        # infinite-distance/no parallax answer
        newrep = UnitSphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
                                             lon=u.Quantity(i_ra, u.radian, copy=False),
                                             copy=False)
    else:
        # When there is a distance, apply the parallax/offset to the SSB as the
        # last step - ensures round-tripping with the icrs_to_cirs transform

        # the distance in intermedrep is *not* a real distance as it does not
        # include the offset back to the SSB
        intermedrep = SphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
                                              lon=u.Quantity(i_ra, u.radian, copy=False),
                                              distance=cirs_coo.distance,
                                              copy=False)

        newxyz = intermedrep.to_cartesian().xyz
        # roll xyz to last axis and add the barycentre position
        newxyz = np.rollaxis(newxyz, 0, newxyz.ndim) + astrom['eb'] * u.au
        # roll xyz back to the first axis
        newxyz = np.rollaxis(newxyz, -1, 0)
        newrep = CartesianRepresentation(newxyz).represent_as(SphericalRepresentation)

    return icrs_frame.realize_frame(newrep)