def merge_cells(angle_indices):
    def connection(index_a, index_b):
        if (index_a == 0 or index_a == 2) and (index_b == 0 or index_b == 2):
            angle_difference = min((index_a - index_b) % 180, (index_b - index_a) % 180)
            return abs(angle_difference) <= 0
        else:
            return False

    grid_size = len(angle_indices)
    num_grid_elements = grid_size * grid_size

    graph = np.zeros((num_grid_elements, num_grid_elements))
    for i in range(grid_size):
        for j in range(grid_size):
            index = np.ravel_multi_index((i, j), (grid_size, grid_size))

            if i > 0 and connection(angle_indices[i, j], angle_indices[i - 1, j]):
                graph[index, np.ravel_multi_index((i - 1, j), (grid_size, grid_size))] = 1

            if i < (grid_size - 1) and connection(angle_indices[i, j], angle_indices[i + 1, j]):
                graph[index, np.ravel_multi_index((i + 1, j), (grid_size, grid_size))] = 1

            if j > 0 and connection(angle_indices[i, j], angle_indices[i, j - 1]):
                graph[index, np.ravel_multi_index((i, j - 1), (grid_size, grid_size))] = 1

            if j < (grid_size - 1) and connection(angle_indices[i, j], angle_indices[i, j + 1]):
                graph[index, np.ravel_multi_index((i, j + 1), (grid_size, grid_size))] = 1

    n_components, labels = connected_components(graph, directed=False, return_labels=True)
    return (labels + 1).reshape(grid_size, grid_size)

def generateGridAdj(nrows, ncols):
    Adj = np.zeros((nrows * ncols,nrows * ncols))
    
    
    for i in range(ncols):
        for j in range(nrows):
            k = np.ravel_multi_index((i,j), dims=(ncols, nrows), order='F')
            if i > 0:
                ii = i-1
                jj = j
                kk = np.ravel_multi_index((ii,jj), dims=(ncols, nrows), order='F')
                Adj[k,kk] = 1
                        
            if i<ncols-1:
                ii=i+1
                jj=j
                kk = np.ravel_multi_index((ii,jj), dims=(ncols, nrows), order='F')
                Adj[k,kk] = 1
            
            if j>0:
                ii=i
                jj=j-1
                kk = np.ravel_multi_index((ii,jj), dims=(ncols, nrows), order='F')
                Adj[k,kk] = 1
            if j<nrows-1:
                ii=i
                jj=j+1
                kk = np.ravel_multi_index((ii,jj), dims=(ncols, nrows), order='F')
                Adj[k,kk] = 1
        
    return Adj

def create_sparse_pre_rate_matrix(state_space_shape, row, col, rate):
    """
    Create the pre-rate matrix with empty diagonal.

    """
    # check conformability of input arrays
    ndim = len(state_space_shape)
    assert_equal(len(row.shape), 2)
    assert_equal(len(col.shape), 2)
    assert_equal(len(rate.shape), 1)
    assert_equal(row.shape[0], rate.shape[0])
    assert_equal(col.shape[0], rate.shape[0])
    assert_equal(row.shape[1], ndim)
    assert_equal(col.shape[1], ndim)

    # Define the transition rate matrix after the multivariate process
    # is collapsed to a univariate process.
    nstates = np.prod(state_space_shape)
    mrow = np.ravel_multi_index(row.T, state_space_shape)
    mcol = np.ravel_multi_index(col.T, state_space_shape)

    # Diagonal entries are allowed for computing dwell times,
    # but they are not allowed for transition expectations.
    #assert_(not np.any(mrow == mcol))

    # create the sparse pre_Q matrix from the sparse arrays
    return coo_matrix((rate, (mrow, mcol)), (nstates, nstates))

def TwoPeriodSun(constant,delta,slope,slopedir,lat):
    # First derive A1 and A2 from the normal procedure    
    A1,A2 = SunHours(delta,slope,slopedir,lat)
    
    # Then calculate the other two functions.
    # Initialize function    
    
    a,b,c = Constants(delta,slope,slopedir,lat)
    riseSlope, setSlope = BoundsSlope(a,b,c)  
    
    B1 = np.maximum(riseSlope,setSlope)
    B2 = np.minimum(riseSlope,setSlope)
    
    Angle_B1 = AngleSlope(a,b,c,B1)
    Angle_B2 = AngleSlope(a,b,c,B2) 
    
    B1[abs(Angle_B1) > 0.001] = np.pi - B1[abs(Angle_B1) > 0.001]
    B2[abs(Angle_B2) > 0.001] = -np.pi - B2[abs(Angle_B2) > 0.001]    
    
    # Check if two periods really exist
    ID = np.ravel_multi_index(np.where(np.logical_and(B2 >= A1, B1 >= A2) == True),a.shape)
    Val = IntegrateSlope(constant,B2.flat[ID],B1.flat[ID],delta,slope.flat[ID],slopedir.flat[ID],lat.flat[ID])
    ID = ID[Val < 0]
    
    # Finally calculate resulting values
    Vals = np.zeros(B1.shape)

    Vals.flat[ID] = (IntegrateSlope(constant,A1.flat[ID],B2.flat[ID],delta,slope.flat[ID],slopedir.flat[ID],lat.flat[ID])  + 
                       IntegrateSlope(constant,B1.flat[ID],A2.flat[ID],delta,slope.flat[ID],slopedir.flat[ID],lat.flat[ID]))
    ID = np.ravel_multi_index(np.where(Vals == 0),a.shape)   
    Vals.flat[ID] = IntegrateSlope(constant,A1.flat[ID],A2.flat[ID],delta,slope.flat[ID],slopedir.flat[ID],lat.flat[ID])
    
    return(Vals)

 def compute_offsets(self):
 
 #    Compute neighbor offsets index.
     mask_coords = self.mask_coords
     vol_num = self.vol_num
     nb_num = self.nb_num
     dims = self.dims
     mask_idx = np.ravel_multi_index(self.mask_coords.T, dims)
      
     offsets = self.nb_shape.compute_off()
     ds_idx = np.arange(vol_num)
     volid = np.zeros(nb_num*vol_num)
     nbid = np.zeros(nb_num*vol_num)
     is_nb = np.zeros(nb_num*vol_num)
     
     
     count = 0
     for v in range(self.vol_num):
         v3d = mask_coords[v,:]
         nb3d = v3d + offsets.T
         imgnb = is_in_image(nb3d, dims)
         nb1d = np.ravel_multi_index(nb3d[imgnb, :].T, dims)
         is_in_mask = np.in1d(nb1d,mask_idx)
         nb1d = nb1d[is_in_mask]
         nb_idx_in_mask = np.in1d(mask_idx,nb1d)
         nb1d = mask_idx[nb_idx_in_mask]
         nb_idx_in_ds = np.nonzero(nb_idx_in_mask)[0]
         volnb_num = len(nb1d)
         volid[ count: count + volnb_num]= np.tile(v, volnb_num)
         nbid[ count: count + volnb_num] = nb_idx_in_ds
         is_nb[ count: count + volnb_num] = 1
         count  = count + volnb_num
     neighbor_sparse_matrix = sparse.csc_matrix((is_nb,(volid,nbid)),shape = (vol_num,vol_num))
     return neighbor_sparse_matrix

def _offsets_to_raveled_neighbors(image_shape, structure, center):
    """Compute offsets to a samples neighbors if the image would be raveled.

    Parameters
    ----------
    image_shape : tuple
        The shape of the image for which the offsets are computed.
    structure : ndarray
        A structuring element determining the neighborhood expressed as an
        n-D array of 1's and 0's.
    center : sequence
        Tuple of indices specifying the center of `selem`.

    Returns
    -------
    offsets : ndarray
        Linear offsets to a samples neighbors in the raveled image, sorted by
        their Euclidean distance from the center.

    Examples
    --------
    >>> _offsets_to_raveled_neighbors((4, 5), np.ones((4, 3)), (1, 1))
    array([-5, -1,  1,  5, -6, -4,  4,  6, 10,  9, 11])
    """
    structure = structure.copy()  # Don't modify original input
    structure[tuple(center)] = 0  # Ignore the center; it's not a neighbor
    connection_indices = np.transpose(np.nonzero(structure))
    offsets = (np.ravel_multi_index(connection_indices.T, image_shape) -
               np.ravel_multi_index(center, image_shape))
    squared_distances = np.sum((connection_indices - center) ** 2, axis=1)
    return offsets[np.argsort(squared_distances)]

        def setInSlot(self, slot, key, value):
            shape = self.inputs["shape"].value
            eraseLabel = self.inputs["eraser"].value
            neutralElement = 0
    
            self.lock.acquire()
            #fix slicing of single dimensions:            
            start, stop = sliceToRoi(key, shape, extendSingleton = False)
            start = start.floor()
            stop = stop.floor()
            
            tempKey = roiToSlice(start-start, stop-start, hardBind = True)
            
            stop += numpy.where(stop-start == 0,1,0)

            key = roiToSlice(start,stop)

            updateShape = tuple(stop-start)
    
            update = self._denseArray[key].copy()
            
            update[tempKey] = value

            startRavel = numpy.ravel_multi_index(numpy.array(start, numpy.int32),shape)
            
            #insert values into dict
            updateNZ = numpy.nonzero(numpy.where(update != neutralElement,1,0))
            updateNZRavelSmall = numpy.ravel_multi_index(updateNZ, updateShape)
            
            if isinstance(value, numpy.ndarray):
                valuesNZ = value.ravel()[updateNZRavelSmall]
            else:
                valuesNZ = value
    
            updateNZRavel = numpy.ravel_multi_index(updateNZ, shape)
            updateNZRavel += startRavel        
    
            self._denseArray.ravel()[updateNZRavel] = valuesNZ        
            
            valuesNZ = self._denseArray.ravel()[updateNZRavel]
            
            self._denseArray.ravel()[updateNZRavel] =  valuesNZ       
    
            
            td = blist.sorteddict(zip(updateNZRavel.tolist(),valuesNZ.tolist()))
       
            self._sparseNZ.update(td)
            
            #remove values to be deleted
            updateNZ = numpy.nonzero(numpy.where(update == eraseLabel,1,0))
            if len(updateNZ)>0:
                updateNZRavel = numpy.ravel_multi_index(updateNZ, shape)
                updateNZRavel += startRavel    
                self._denseArray.ravel()[updateNZRavel] = neutralElement
                for index in updateNZRavel:
                    self._sparseNZ.pop(index)
            
            self.lock.release()
            
            self.outputs["Output"].setDirty(key)

def reverse_interpolate_two_array(value1, array1, value2, array2, delta1=0.1, delta2=0.1):
    """
    Tries to reverse interpolate two vales from two arrays with the same dimensions, and finds a common index
    for value1 and value2 in their respective arrays. the deltas define the search radius for a close value match
    to the arrays.

    :return: index1, index2
    """
    tth_ind = np.argwhere(np.abs(array1 - value1) < delta1)
    azi_ind = np.argwhere(np.abs(array2 - value2) < delta2)

    tth_ind_ravel = np.ravel_multi_index((tth_ind[:, 0], tth_ind[:, 1]), dims=array1.shape)
    azi_ind_ravel = np.ravel_multi_index((azi_ind[:, 0], azi_ind[:, 1]), dims=array2.shape)

    common_ind_ravel = np.intersect1d(tth_ind_ravel, azi_ind_ravel)
    result_ind = np.unravel_index(common_ind_ravel, dims=array1.shape)

    while len(result_ind[0]) > 1:
        if np.max(np.diff(array1)) > 0:
            delta1 = np.max(np.diff(array1[result_ind]))

        if np.max(np.diff(array2)) > 0:
            delta2 = np.max(np.diff(array2[result_ind]))

        tth_ind = np.argwhere(np.abs(array1[result_ind] - value1) < delta1)
        azi_ind = np.argwhere(np.abs(array2[result_ind] - value2) < delta2)

        print(result_ind)

        common_ind = np.intersect1d(tth_ind, azi_ind)
        result_ind = (result_ind[0][common_ind], result_ind[1][common_ind])

    return result_ind[0], result_ind[1]

def test_rmi():
    I1 = _rmi([3, 4], [10, 10])
    assert_equal(I1, 34)
    I1 = _rmi([0, 0], [10, 10])
    assert_equal(I1, 0)
    assert_raises(ValueError, _rmi, [10, 0], [10, 10])

    try:
        from numpy import ravel_multi_index
    except ImportError:
        raise nose.SkipTest()

    # Dtype of random integers is system dependent
    A, B, C, D = np.random.randint(0, 1000, size=[4, 100])
    I1 = _rmi([A, B], dims=[1000, 1000])
    I2 = ravel_multi_index([A, B], dims=[1000, 1000])
    assert_array_equal(I1, I2)
    I1 = _rmi([A, B, C, D], dims=[1000]*4)
    I2 = ravel_multi_index([A, B, C, D], dims=[1000]*4)
    assert_array_equal(I1, I2)
    # Check for overflow with small int types
    indices = np.random.randint(0, 255, size=(2, 100))
    dims = (1000, 1000)
    I1 = _rmi(indices, dims=dims)
    I2 = ravel_multi_index(indices, dims=dims)
    assert_array_equal(I1, I2)

    def __init__(self, model):
        """ Initialize the filtered stim model
        """
        self.model = model
        N = model['N']
        prms = model['network']['weight']

        self.prior = create_prior(prms['prior'])
	
        # Implement refractory period by having negative mean on self loops
        if 'refractory_prior' in prms:
            #self.mu[np.diag_indices(N)] = prms['mu_refractory']
            self.refractory_prior = create_prior(prms['refractory_prior'])

            # Get the upper and lower diagonal indices so that we can evaluate
            # the log prob of the refractory weights separately from the
            # log prob of the regular weights
            self.diags = np.ravel_multi_index(np.diag_indices(N), (N,N))
            lower = np.ravel_multi_index(np.tril_indices(N,k=-1), (N,N))
            upper = np.ravel_multi_index(np.triu_indices(N,k=1), (N,N))
            self.nondiags = np.concatenate((lower, upper))

        # Define weight matrix
        self.W_flat = T.dvector(name='W')
        self.W = T.reshape(self.W_flat,(N,N))

        if hasattr(self, 'refractory_prior'):
            self.log_p = self.prior.log_p(self.W.take(self.nondiags)) + \
                         self.refractory_prior.log_p(self.W.take(self.diags))
        else:
            self.log_p = self.prior.log_p(self.W)

    def _init_indices(self):
        # Moore neighbourhood for now. Will see others later.
        self.n_neighbours = 8 # Moore
        cols, rows = map(np.ndarray.flatten, np.meshgrid(*self._valid_range()))
        self.inner_indices = np.ravel_multi_index((rows, cols),
                                                  self.dims,
                                                  order='C')
        left_cols, right_cols = cols - 1, cols + 1
        up_rows, down_rows = rows - 1, rows + 1
        neighbour_pix_locs = [(up_rows, left_cols), (up_rows, cols),
                              (up_rows, right_cols),
                              (rows, left_cols), (rows, cols),
                              (rows, right_cols),
                              (down_rows, left_cols), (down_rows, cols),
                              (down_rows, right_cols)]

        self.indices = np.empty((self.inner_indices.size,
                                 self.n_neighbours + 1),
                                dtype=np.int)
        for i in xrange(len(neighbour_pix_locs)):
            self.indices[:,i] \
                    = np.ravel_multi_index(neighbour_pix_locs[i],
                                           self.dims,
                                           mode='wrap',
                                           order='C')

    def test_dtypes(self):
        # Test with different data types
        for dtype in [np.int16, np.uint16, np.int32,
                      np.uint32, np.int64, np.uint64]:
            coords = np.array(
                [[1, 0, 1, 2, 3, 4], [1, 6, 1, 3, 2, 0]], dtype=dtype)
            shape = (5, 8)
            uncoords = 8*coords[0]+coords[1]
            assert_equal(np.ravel_multi_index(coords, shape), uncoords)
            assert_equal(coords, np.unravel_index(uncoords, shape))
            uncoords = coords[0]+5*coords[1]
            assert_equal(
                np.ravel_multi_index(coords, shape, order='F'), uncoords)
            assert_equal(coords, np.unravel_index(uncoords, shape, order='F'))

            coords = np.array(
                [[1, 0, 1, 2, 3, 4], [1, 6, 1, 3, 2, 0], [1, 3, 1, 0, 9, 5]],
                dtype=dtype)
            shape = (5, 8, 10)
            uncoords = 10*(8*coords[0]+coords[1])+coords[2]
            assert_equal(np.ravel_multi_index(coords, shape), uncoords)
            assert_equal(coords, np.unravel_index(uncoords, shape))
            uncoords = coords[0]+5*(coords[1]+8*coords[2])
            assert_equal(
                np.ravel_multi_index(coords, shape, order='F'), uncoords)
            assert_equal(coords, np.unravel_index(uncoords, shape, order='F'))

def create_sparse_rate_matrix(state_space_shape, row, col, rate):
    """
    Create the rate matrix.

    """
    # check conformability of input arrays
    ndim = len(state_space_shape)
    assert_equal(len(row.shape), 2)
    assert_equal(len(col.shape), 2)
    assert_equal(len(rate.shape), 1)
    assert_equal(row.shape[0], rate.shape[0])
    assert_equal(col.shape[0], rate.shape[0])
    assert_equal(row.shape[1], ndim)
    assert_equal(col.shape[1], ndim)

    # create the sparse Q matrix from the sparse arrays
    nstates = np.prod(state_space_shape)
    mrow = np.ravel_multi_index(row.T, state_space_shape)
    mcol = np.ravel_multi_index(col.T, state_space_shape)
    Q = coo_matrix((rate, (mrow, mcol)), (nstates, nstates))

    # get the dense array of exit rates, and set the diagonal
    exit_rates = Q.sum(axis=1).A.flatten()
    Q.setdiag(-exit_rates)

    return Q

def dig(this, next):
    # first select active mesh from next level, aligned to this level
    next.data = next.data[activemask(next, this.Nmesh)]

    # calculate the index of the covered mesh on this level
    index_from_next = numpy.ravel_multi_index(
        next.data["gps"].T // (next.Nmesh // this.Nmesh), (this.Nmesh, this.Nmesh, this.Nmesh)
    )

    # index this level
    index = numpy.ravel_multi_index(this.data["gps"].T, (this.Nmesh, this.Nmesh, this.Nmesh))

    preserve_mask2 = exist(index, index_from_next)
    next.data = next.data[preserve_mask2]

    if next.DownSample > 1:
        # downsample levels need to carry on the interpolation
        # of previous level
        next.data["delta"] += trilinearinterp(
            next.data["gps"] * (1.0 * this.Nmesh / next.Nmesh),
            this.data["gps"],
            this.data["delta"],
            (next.Nmesh, next.Nmesh, next.Nmesh),
        )
        next.data["disp"] += trilinearinterp(
            next.data["gps"] * (1.0 * this.Nmesh / next.Nmesh),
            this.data["gps"],
            this.data["disp"],
            (next.Nmesh, next.Nmesh, next.Nmesh),
        )

    # locate those to be removed
    covered_mask = exist(index_from_next, index)
    preserve_mask = ~covered_mask
    this.data = this.data[preserve_mask]

 def map_to_phon(self):
     '''decibel are converted into phon which are a logarithmic unit to human loudness sensation'''
     # number of bark bands, matrix length in time dim
     n_bands = self.processed.shape[0]
     t       = self.processed.shape[1]
 
     
     # DB-TO-PHON BARK-SCALE-LIMIT TABLE
     # introducing 1 level more with level(1) being infinite
     # to avoid (levels - 1) producing errors like division by 0
     table_dim = n_bands; 
     cbv       = np.concatenate((np.tile(np.inf,(table_dim,1)), 
                                 self.loudn_bark[:,0:n_bands].transpose()),1)
     
     # the matrix 'levels' stores the correct Phon level for each datapoint
     # init lowest level = 2
     levels = np.tile(2,(n_bands,t)) 
     
     for lev in range(1,6): 
         db_thislev = np.tile(np.asarray([cbv[:,lev]]).transpose(),(1,t))
         levels[np.where(self.processed > db_thislev)] = lev + 2
     #ok here we compute indexes that match the cbv array such that when indexing the cbv array we get back an matrix with the proper dimensions
     #this indices match the upper or lower phon level according to the current value of spl in our matrix
     cbv_ind_hi = np.ravel_multi_index(dims=(table_dim,7), multi_index=np.array([np.tile(np.array([range(0,table_dim)]).transpose(),(1,t)), levels-1]), order='F') 
     cbv_ind_lo = np.ravel_multi_index(dims=(table_dim,7), multi_index=np.array([np.tile(np.array([range(0,table_dim)]).transpose(),(1,t)), levels-2]), order='F') 
 
     #for each datapoint in our matrix a interpolation factor 0 < f < 1 is calculated
     ifac = (self.processed[:,0:t] - cbv.transpose().ravel()[cbv_ind_lo]) / (cbv.transpose().ravel()[cbv_ind_hi] - cbv.transpose().ravel()[cbv_ind_lo])
     
     ifac[np.where(levels==2)] = 1 # keeps the upper phon value;
     ifac[np.where(levels==8)] = 1 # keeps the upper phon value;
 
     #finally the interpolation is computed here
     self.processed[:,0:t] = AudioAnalytics.phons.transpose().ravel()[levels - 2] + (ifac * (AudioAnalytics.phons.transpose().ravel()[levels - 1] - AudioAnalytics.phons.transpose().ravel()[levels - 2])) # OPT: pre-calc diff    

def transform2phon(matrix):

    # number of bark bands, matrix length in time dim
    n_bands = matrix.shape[0]
    t       = matrix.shape[1]

    # DB-TO-PHON BARK-SCALE-LIMIT TABLE
    # introducing 1 level more with level(1) being infinite
    # to avoid (levels - 1) producing errors like division by 0

    #%%table_dim = size(CONST_loudn_bark,2);
    table_dim = n_bands; # OK
    cbv       = np.concatenate((np.tile(np.inf,(table_dim,1)), loudn_bark[:,0:n_bands].transpose()),1) # OK

    # init lowest level = 2
    levels = np.tile(2,(n_bands,t)) # OK

    for lev in range(1,6): # OK
        db_thislev = np.tile(np.asarray([cbv[:,lev]]).transpose(),(1,t))
        levels[np.where(matrix > db_thislev)] = lev + 2

    # the matrix 'levels' stores the correct Phon level for each data point
    cbv_ind_hi = np.ravel_multi_index(dims=(table_dim,7), multi_index=np.array([np.tile(np.array([range(0,table_dim)]).transpose(),(1,t)), levels-1]), order='F')
    cbv_ind_lo = np.ravel_multi_index(dims=(table_dim,7), multi_index=np.array([np.tile(np.array([range(0,table_dim)]).transpose(),(1,t)), levels-2]), order='F')

    # interpolation factor % OPT: pre-calc diff
    ifac = (matrix[:,0:t] - cbv.transpose().ravel()[cbv_ind_lo]) / (cbv.transpose().ravel()[cbv_ind_hi] - cbv.transpose().ravel()[cbv_ind_lo])

    ifac[np.where(levels==2)] = 1 # keeps the upper phon value;
    ifac[np.where(levels==8)] = 1 # keeps the upper phon value;

    # phons has been initialized globally above

    matrix[:,0:t] = phons.transpose().ravel()[levels - 2] + (ifac * (phons.transpose().ravel()[levels - 1] - phons.transpose().ravel()[levels - 2])) # OPT: pre-calc diff
    return(matrix)

def unpack_array_c(input_data, output_shape):
    """Unpack an array."""
    input_data  = np.asarray(input_data)
    output_data = _prepare_arrays(input_data, output_shape)
    input_shape = input_data.shape
    
    # Total number of elements in each type of array.
    nft = np.prod(input_shape)
    nn  = np.prod(output_shape)
    
    # Extract single dimension parameters, where relevant.
    ny = output_shape[0]
    nf = input_shape[-1]
    nx = output_shape[-1]
    
    # Determine the central element.
    nfo = nf - 1 if nx % 2 == 0 else nf
    for i in range(nft):
        
        # Extract the cartesian coordinates in the original grid.
        y = i // nf
        x = i %  nf
        
        # re-map onto the wider grid.
        io = (y * nx) + x
        ii = (y * nf) + x
        output_data.flat[io] = input_data.flat[ii]
        
        # Check our indexing math with numpy's version.
        assert np.allclose([y,x], np.unravel_index(ii, input_shape))
        assert np.allclose([io], np.ravel_multi_index([y, x], output_shape))
        
        # If we are at y=0, flip around center.
        if y == 0 and 0 < x < nfo:
            yo = 0
            xo = nx - x
            io = (yo * nx) + xo
            
            # Check our math!
            assert np.allclose([io], np.ravel_multi_index([yo, xo], output_shape))
            assert np.allclose([yo, xo], np.unravel_index(io, output_shape))
            
            # Insert data.
            output_data.flat[io] = np.conj(input_data.flat[ii])
            
        # If we are beyond y=0, flip in both axes.
        if y > 0 and 0 < x < nfo:
            yo = ny - y
            xo = nx - x
            io = (yo * nx) + xo
            
            # Check our math!
            assert np.allclose([io], np.ravel_multi_index([yo, xo], output_shape))
            assert np.allclose([yo, xo], np.unravel_index(io, output_shape))
            
            # Insert data.
            output_data.flat[io] = np.conj(input_data.flat[ii])
    
    return output_data

 def test_clipmodes(self):
     # Test clipmodes
     assert_equal(np.ravel_multi_index([5,1,-1,2], (4,3,7,12), mode='wrap'),
                  np.ravel_multi_index([1,1,6,2], (4,3,7,12)))
     assert_equal(np.ravel_multi_index([5,1,-1,2], (4,3,7,12),
                              mode=('wrap','raise','clip','raise')),
                  np.ravel_multi_index([1,1,0,2], (4,3,7,12)))
     assert_raises(ValueError, np.ravel_multi_index, [5,1,-1,2], (4,3,7,12))

def build_diff_matrix(int_band, dxvec, shape):
    """Builds a Laplacian matrix"""
    # TODO: generalize this!
    shape = shape.astype(np.int64)
    dim = len(shape)
    if dim == 3:
        # Assume now that band is not a linear index, but a 3D index (ie, int_band)
        dx, dy, dz = dxvec
        weights = np.array([-2/dx**2 -2/dy**2 -2/dz**2,
                             1/dx**2, 1/dx**2,
                             1/dy**2, 1/dy**2,
                             1/dz**2, 1/dz**2])

        offsets = np.array([[ 0, 0, 0],
                            [ 1, 0, 0],
                            [-1, 0, 0],
                            [ 0, 1, 0],
                            [ 0,-1, 0],
                            [ 0, 0, 1],
                            [ 0, 0,-1]])

        Li = np.tile(np.arange(int_band.shape[0])[:, np.newaxis], weights.size)
        Lj = np.zeros_like(Li)
        Ls = np.zeros_like(Li, dtype=np.float)
        i, j, k = int_band.T  #np.ravel_multi_index(band.T, shape)

        for c in xrange(weights.size):
            ii = i + offsets[c, 0]
            jj = j + offsets[c, 1]
            kk = k + offsets[c, 2]
            Ls[:, c] = weights[c]
            Lj[:, c] = np.ravel_multi_index((ii, jj, kk), shape)
        L = coo_matrix((Ls.ravel(), (Li.ravel(), Lj.ravel())),
                       (int_band.shape[0], shape[0] * shape[1] * shape[2])).tocsr()
        return L[:, np.ravel_multi_index(int_band.T, shape)]
    elif dim == 2:
        dx, dy, dz = dxvec
        weights = np. array([-2/dx**2 -2/dy**2,
                              1/dx**2, 1/dx**2,
                              1/dy**2, 1/dy**2])
        offsets = np.array([[ 0, 0],
                            [ 1, 0],
                            [-1, 0],
                            [ 0, 1],
                            [ 0,-1]])
        Li = np.tile(np.arange(int_band.shape[0])[:, np.newaxis], weights.size)
        Lj = np.zeros_like(Li)
        Ls = np.zeros_like(Li, dtype=np.float)
        i, j = int_band.T  #np.ravel_multi_index(band.T, shape)

        for c in xrange(weights.size):
            ii = i + offsets[c, 0]
            jj = j + offsets[c, 1]
            Ls[:, c] = weights[c]
            Lj[:, c] = np.ravel_multi_index((ii, jj), shape)
        L = coo_matrix((Ls.ravel(), (Li.ravel(), Lj.ravel())),
                       (int_band.shape[0], shape[0] * shape[1])).tocsr()
        return L[:, np.ravel_multi_index(int_band.T, shape)]

    def _make_hdu_sparse(data, npix, hdu, header):
        shape = data.shape

        # We make a copy, because below we modify `data` to handle non-finite entries
        # TODO: The code below could probably be simplified to use expressions
        # that create new arrays instead of in-place modifications
        # But first: do we want / need the non-finite entry handling at all and always cast to 64-bit float?
        data = data.copy()

        if len(shape) == 2:
            data_flat = np.ravel(data)
            data_flat[~np.isfinite(data_flat)] = 0
            nonzero = np.where(data_flat > 0)
            value = data_flat[nonzero].astype(float)

            cols = [
                fits.Column('PIX', 'J', array=nonzero[0]),
                fits.Column('VALUE', 'E', array=value),
            ]
        elif npix[0].size == 1:
            shape_flat = shape[:-2] + (shape[-1] * shape[-2],)
            data_flat = np.ravel(data).reshape(shape_flat)
            data_flat[~np.isfinite(data_flat)] = 0
            nonzero = np.where(data_flat > 0)
            channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])
            value = data_flat[nonzero].astype(float)

            cols = [
                fits.Column('PIX', 'J', array=nonzero[-1]),
                fits.Column('CHANNEL', 'I', array=channel),
                fits.Column('VALUE', 'E', array=value),
            ]
        else:
            data_flat = []
            channel = []
            pix = []
            for i, _ in np.ndenumerate(npix[0]):
                data_i = np.ravel(data[i[::-1]])
                data_i[~np.isfinite(data_i)] = 0
                pix_i = np.where(data_i > 0)
                data_i = data_i[pix_i]
                data_flat += [data_i]
                pix += pix_i
                channel += [np.ones(data_i.size, dtype=int) *
                            np.ravel_multi_index(i[::-1], shape[:-2])]

            pix = np.concatenate(pix)
            channel = np.concatenate(channel)
            value = np.concatenate(data_flat).astype(float)

            cols = [
                fits.Column('PIX', 'J', array=pix),
                fits.Column('CHANNEL', 'I', array=channel),
                fits.Column('VALUE', 'E', array=value),
            ]

        return fits.BinTableHDU.from_columns(cols, header=header, name=hdu)