def sparseScalarProductOfDot (A, B, C, out=None):
    '''
    Returns A * np.dot(B, C), however it does so keeping in  mind 
    the sparsity of A, calculating values only where required.
     
    Params
    A         - a sparse CSR matrix
    B         - a dense matrix
    C         - a dense matrix
    out       - if specified, must be a sparse CSR matrix with identical
                non-zero pattern to A (i.e. same indices and indptr)
    
    Returns
    out_data, though note that this is the same parameter passed in and overwitten.
    '''
    assert ssp.isspmatrix_csr(A), "A matrix is not a CSR matrix"
    assert not np.isfortran(B), "B matrix is not stored in row-major form"
    assert not np.isfortran(C), "C matrix is not stored in row-major form"

    if out is None:
        out = A.copy()
    if A.dtype == np.float64:
        compiled.sparseScalarProductOfDot_f8(A.data, A.indices, A.indptr, B, C, out.data)
    elif A.dtype == np.float32:
        compiled.sparseScalarProductOfDot_f4(A.data, A.indices, A.indptr, B, C, out.data)
    else:
        _sparseScalarProductOfDot_py(A,B,C, out)
    return out

def generate_itp_pn(lut):
    ndim=lut.ndim-1
    
    if ndim == 1: interpolator = interpolators.interpol_1pn
    elif ndim == 2: interpolator = interpolators.interpol_2pn
    elif ndim == 3: interpolator = interpolators.interpol_3pn
    elif ndim == 4: interpolator = interpolators.interpol_4pn
    elif ndim == 5: interpolator = interpolators.interpol_5pn
    elif ndim == 6: interpolator = interpolators.interpol_6pn
    elif ndim >= 7: interpolator = interpolators.interpol_npn
    else:
        print 'Not implemented'
        return
    
    if ndim <= 6:
        if np.isfortran(lut): my_lut=lut
        else: my_lut=np.asfortranarray(lut)        
        def function(wo): 
            return interpolator(wo,my_lut)
    else: 
        if np.isfortran(lut): 
            flat_lut=np.asfortranarray(lut.reshape((-1,lut.shape[-1]) ,order='C'))
        else: 
            flat_lut=np.asfortranarray(lut.reshape((-1,lut.shape[-1]), order='C'))
        shape=np.array(lut.shape)
        size=shape[:-1].prod()
        def function(wo):
            return interpolator(wo,flat_lut,shape[0:-1],shape[-1])
    
    return function

def generate_nan_itp(lut):
    ndim=lut.ndim
    
    if ndim == 1: interpolator = interpolators.interpol_nan_1
    elif ndim == 2: interpolator = interpolators.interpol_nan_2
    elif ndim == 3: interpolator = interpolators.interpol_nan_3
    elif ndim == 4: interpolator = interpolators.interpol_nan_4
    elif ndim == 5: interpolator = interpolators.interpol_5
    elif ndim == 6: interpolator = interpolators.interpol_6
    elif ndim == 7: interpolator = interpolators.interpol_7
    elif ndim >= 8: interpolator = interpolators.interpol_n
    else:
        print 'Not implemented'
        return
    
    if ndim <= 7:
        if np.isfortran(lut): my_lut=lut.astype(float)
        else: my_lut=np.asfortranarray(lut.astype(float))        
        def function(wo): 
            return interpolator(wo,my_lut)
    else: 
        if np.isfortran(lut): flat_lut=lut.astype(float).ravel('C')
        else: flat_lut=lut.astype(float).T.ravel('F')
        shape=np.array(lut.shape)
        size=lut.size
        def function(wo):
            return interpolator(wo,flat_lut,shape)
    
    return function

def remove_adjacencies(labels, bwconn):
    test = np.zeros((2,2), dtype = np.uint32)
    if type(labels) != type(test):
        raise Exception('In remove_adjacencies, labels is not a *NumPy* array')
    if len(labels.shape) != 3:
        raise Exception('In remove_adjacencies, labels is not 3 dimensional')
    if not labels.flags.contiguous or np.isfortran(labels):
        raise Exception('In remove_adjacencies, labels not contiguous or not C-order')
    if labels.dtype != test.dtype:
        raise Exception('In remove_adjacencies, labels not uint32')

    testB=np.zeros((2,2),dtype=np.bool)
    if type(bwconn) != type(testB):
        raise Exception( 'In remove_adjacencies, bwconn is not *NumPy* array')
    if len(bwconn.shape) != 3:
        raise Exception( 'In remove_adjacencies, bwconn is not 3 dimensional')
    if not bwconn.flags.contiguous or np.isfortran(bwconn):
        raise Exception( 'In remove_adjacencies, bwconn not C-order contiguous')
    if bwconn.dtype != testB.dtype:
        raise Exception( 'In remove_adjacencies, bwconn not correct data type')

    sz =  [x+2 for x in labels.shape]   # need border for neighborhoods around the edge voxels
    lbls = np.zeros(sz, dtype=labels.dtype); lbls[1:-1,1:-1,1:-1] = labels
    _pyCext.remove_adjacencies(lbls, bwconn)
    return lbls[1:-1,1:-1,1:-1]

def label_affinities(affinities, labels, nextlabel, threshold):
    test=np.zeros((2,2)); testI=np.zeros((2,2),dtype=np.uint32)
    if type(affinities) != type(test):
        raise Exception( 'In label_affinities, affinities is not *NumPy* array')
    if len(affinities.shape) != 4:
        raise Exception( 'In label_affinities, affinities shape not 4 dimensional')
    if affinities.shape[3] != 3:
        raise Exception( 'In label_affinities, affinities not 3 dimensional')
    if not affinities.flags.contiguous or np.isfortran(affinities):
        raise Exception( 'In label_affinities, affinities not C-order contiguous')
    if affinities.dtype != np.single:
        raise Exception( 'In label_affinities, affinities not single floats')
    if type(labels) != type(testI):
        raise Exception( 'In label_affinities, labels is not *NumPy* array')
    if len(labels.shape) != 3:
        raise Exception( 'In label_affinities, labels is not 3 dimensional')
    if not labels.flags.contiguous or np.isfortran(labels):
        raise Exception( 'In label_affinities, labels not C-order contiguous')
    if labels.dtype != np.uint32:
        raise Exception( 'In label_affinities, labels not uint32')
    if type(threshold) != type(1.0):
        raise Exception( 'In label_affinities, threshold argument is not a float')
    if type(nextlabel) != type(1):
        raise Exception( 'In label_affinities, nextlabel argument is not a integer')
    if nextlabel < 1:
        raise Exception( 'In label_components, nextlabel argument is less than one')

    return _pyCext.label_affinities(affinities,labels,nextlabel,threshold)

 def __init__(self, data, block_length=1, use_blocks=None, offsets=None):
     """
     data can be a numpy array (in C order), a tuple of such arrays, or a list of such tuples or arrays
     """
     self.files = list()
     if isinstance(data, list): #Several files
         for file in data:
             if isinstance(file, tuple):
                 for d in file:
                     assert(isinstance(d, np.ndarray) and not np.isfortran(d))
                 self.files.append(file)
     elif isinstance(data, tuple): #Just one file
         for d in data:
             assert(isinstance(d, np.ndarray) and d.ndim == 2 and not np.isfortran(d))
         self.files.append(data)
     elif isinstance(data, np.ndarray): #One file with one kind of element only (not input-output)
         assert(isinstance(data, np.ndarray) and not np.isfortran(data))
         self.files.append(tuple([data]))
     # Support for block datapoints
     self.block_length = block_length
     if block_length == 1:
         self.block_lengths = [np.int(1)] * self.get_arity()
         self.offsets = [np.int(0)] * self.get_arity()
     elif block_length > 1:
         self.block_lengths = [np.int(block_length) if ub else np.int(1) for ub in use_blocks]  # np.asarray(dtype=np.int) and [np.int(x)] have elements with diff type. Careful!
         self.offsets = [np.int(off) for off in offsets]
         for ub, off in zip(use_blocks, offsets):
             if off != 0 and ub:
                 raise Exception("Can't have both a block size greater than 1 and an offset.")
     else:
         raise Exception("Block size must be positive")

def sparseScalarQuotientOfNormedDot (A, B, C, d, out=None):
    '''
    Returns A / np.dot(B, C/D), however it does so keeping in  mind
    the sparsity of A, calculating values only where required.

    Params
    A         - a sparse CSR matrix
    B         - a dense matrix
    C         - a dense matrix
    d         - a dense vector whose dimensionality matches the column-count of C
    out       - if specified, must be a sparse CSR matrix with identical
                non-zero pattern to A (i.e. same indices and indptr)

    Returns
    out_data, though note that this is the same parameter passed in and overwitten.
    '''
    assert ssp.isspmatrix_csr(A), "A matrix is not a CSR matrix"
    assert not np.isfortran(B), "B matrix is not stored in row-major form"
    assert not np.isfortran(C), "C matrix is not stored in row-major form"

    if out is None:
        out = A.copy()

    if A.dtype == np.float64:
        compiled.sparseScalarQuotientOfNormedDot_f8(A.data, A.indices, A.indptr, B, C, d, out.data)
    elif A.dtype == np.float32:
        compiled.sparseScalarQuotientOfNormedDot_f4(A.data, A.indices, A.indptr, B, C, d, out.data)
    else:
        raise ValueError ("No implementation for the datatype " + str(A.dtype))
    return out

 def isfortran(self, array):
   """
     Same as np.isfortran, but works with lists and tuples.
   """
   if isinstance(array,(list,tuple)):
     return all([np.isfortran(val) for val in array])
   else:
     return np.isfortran(array)

def add_dot(c, a, b):
    if use_blas and isinstance(c, numpy.ndarray):
        if numpy.isfortran(c.T):
            scipy.linalg.blas.dgemm(1., b.T, a.T, 1., c.T, overwrite_c=True)
            return c
        elif numpy.isfortran(c):
            scipy.linalg.blas.dgemm(1., a, b, 1., c, overwrite_c=True)
            return c

    c += numpy.dot(a, b)
    return c

def niiToArray(niiImg, c_contiguos=True, canonical=False):
    if canonical:
        rough_data = nib.as_closest_canonical(niiImg).get_data()
    else:
        rough_data = niiImg.get_data()
    print np.isfortran(rough_data)
    #if c_contiguos and np.isfortran(rough_data):
    #    rough_data = rough_data.T
    if c_contiguos and not rough_data.flags['C_CONTIGUOUS']:
        rough_data = rough_data.T
    return rough_data

def test_ndim_indexing(aggregate_all, ndim, order, outsize=10):
    nindices = int(outsize ** ndim)
    outshape = tuple([outsize] * ndim)
    group_idx = np.random.randint(0, outsize, size=(ndim, nindices))
    a = np.random.random(group_idx.shape[1])
    res = aggregate_all(group_idx, a, size=outshape, order=order)
    if ndim > 1 and order == 'F':
        # 1d arrays always return False here
        assert np.isfortran(res)
    else:
        assert not np.isfortran(res)
    assert res.shape == outshape

 def can_add(self, data):
     """Returns True iff data can be stored.
     This usually means it is of the same kind as previously stored arrays.
     """
     # Obviously we could just store arbitrary arrays in some implementations (e.g. NPY)
     # But lets keep jagged contracts...
     template = self.template()
     if template is None:
         return True
     return (template.dtype >= data.dtype and
             data.shape[-1] == template.shape[-1] and
             np.isfortran(data) == np.isfortran(data))

def add_AB_to_C(A, B, C):
    """
    Compute C += AB in-place.

    This uses gemm from whatever BLAS is available.
    MKL requires Fortran ordered arrays to avoid copies.
    Hence we work with transpositions of default c-style arrays.

    This function throws error if computation is not in-place.
    """
    gemm = sl.get_blas_funcs("gemm", (A, B, C))
    assert np.isfortran(C.T) and np.isfortran(A.T) and np.isfortran(B.T)
    D = gemm(1.0, B.T, A.T, beta=1, c=C.T, overwrite_c=1)
    assert D.base is C or D.base is C.base

    def __init__(self, cl_ctx, h0, eta0, u0, v0):
        #Make sure that the data is single precision floating point
        if (not np.issubdtype(h0.dtype, np.float32) or np.isfortran(h0)):
            print "Converting H0"
            h0 = h0.astype(np.float32, order='C')
            
        if (not np.issubdtype(eta0.dtype, np.float32) or np.isfortran(eta0)):
            print "Converting Eta0"
            eta0 = eta0.astype(np.float32, order='C')
            
        if (not np.issubdtype(u0.dtype, np.float32) or np.isfortran(u0)):
            print "Converting U0"
            u0 = u0.astype(np.float32, order='C')
            
        if (not np.issubdtype(v0.dtype, np.float32) or np.isfortran(v0)):
            print "Converting V0"
            v0 = v0.astype(np.float32, order='C')
        
        self.ny, self.nx = h0.shape
        self.nx = self.nx - 2 # Ghost cells
        self.ny = self.ny - 2

        assert(h0.shape == (self.ny+2, self.nx+2))
        assert(eta0.shape == (self.ny+2, self.nx+2))
        assert(u0.shape == (self.ny+2, self.nx+1))
        assert(v0.shape == (self.ny+1, self.nx+2))

        #Upload data to the device
        mf = cl.mem_flags
        self.h0 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=h0)
        
        self.eta0 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=eta0)
        self.eta1 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=eta0)
        
        self.u0 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=u0)
        self.u1 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=u0)
        
        self.v0 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=v0)
        self.v1 = cl.Buffer(cl_ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf=v0)
        
        self.h0_pitch = np.int32(h0.shape[1]*4)
        
        self.eta0_pitch = np.int32(eta0.shape[1]*4)
        self.eta1_pitch = np.int32(eta0.shape[1]*4)
        
        self.u0_pitch = np.int32(u0.shape[1]*4)
        self.u1_pitch = np.int32(u0.shape[1]*4)
        
        self.v0_pitch = np.int32(v0.shape[1]*4)
        self.v1_pitch = np.int32(v0.shape[1]*4)

 def test_fstecr_fstluk_order(self):
     rmn.fstopt(rmn.FSTOP_MSGLVL,rmn.FSTOPI_MSG_CATAST)
     fname = '__rpnstd__testfile2__.fst'
     try:
         os.unlink(fname)
     except:
         pass
     funit = rmn.fstopenall(fname,rmn.FST_RW)
     (ig1,ig2,ig3,ig4) = rmn.cxgaig(self.grtyp,self.xg14[0],self.xg14[1],self.xg14[2],self.xg14[3])
     (ni,nj) = (90,45)
     la = rmn.FST_RDE_META_DEFAULT.copy()
     la.update(
         {'nomvar' : 'LA',
          'typvar' : 'C',
          'ni' : ni,
          'nj' : nj,
          'nk' : 1,
          'grtyp' : self.grtyp,
          'ig1' : ig1,
          'ig2' : ig2,
          'ig3' : ig3,
          'ig4' : ig4
          }
         )
     lo = la.copy()
     lo['nomvar'] = 'LO'
     #Note: For the order to be ok in the FSTD file, order='FORTRAN' is mandatory
     la['d'] = np.empty((ni,nj),dtype=np.float32,order='FORTRAN')
     lo['d'] = np.empty((ni,nj),dtype=np.float32,order='FORTRAN')
     for j in xrange(nj):
         for i in xrange(ni):
             lo['d'][i,j] = 100.+float(i)        
             la['d'][i,j] = float(j)
     rmn.fstecr(funit,la['d'],la)
     rmn.fstecr(funit,lo)
     rmn.fstcloseall(funit)
     funit = rmn.fstopenall(fname,rmn.FST_RW)
     kla = rmn.fstinf(funit,nomvar='LA')['key']
     la2 = rmn.fstluk(kla)#,rank=2)
     klo = rmn.fstinf(funit,nomvar='LO')['key']
     lo2 = rmn.fstluk(klo)#,rank=2)
     rmn.fstcloseall(funit)
     try:
         os.unlink(fname)
     except:
         pass
     self.assertTrue(np.isfortran(la2['d']))
     self.assertTrue(np.isfortran(lo2['d']))
     self.assertFalse(np.any(np.fabs(la2['d'] - la['d']) > self.epsilon))
     self.assertFalse(np.any(np.fabs(lo2['d'] - lo['d']) > self.epsilon))

def test_cmp_ndim(aggregate_cmp, ndim, order, outsize=100, decimal=14):
    nindices = int(outsize ** ndim)
    outshape = tuple([outsize] * ndim)
    group_idx = np.random.randint(0, outsize, size=(ndim, nindices))
    a = np.random.random(group_idx.shape[1])

    res = aggregate_cmp.func(group_idx, a, size=outshape, order=order)
    ref = aggregate_cmp.func_ref(group_idx, a, size=outshape, order=order)
    if ndim > 1 and order == 'F':
        # 1d arrays always return False here
        assert np.isfortran(res)
    else:
        assert not np.isfortran(res)
    assert res.shape == outshape
    np.testing.assert_array_almost_equal(res, ref, decimal=decimal)

def add_outer(a, x, y):
    """Add the outer product of x and y to a, possibly overwriting a.
    a = add_outer(a, x, y) is equivalent to a += numpy.outer(x, y)."""

    if use_blas and isinstance(a, numpy.ndarray):
        if numpy.isfortran(a.T):
            scipy.linalg.blas.dger(1., y, x, a=a.T, overwrite_a=1)
            return a
        elif numpy.isfortran(a):
            scipy.linalg.blas.dger(1., x, y, a=a, overwrite_a=1)
            return a

    # einsum is written in C and is faster than outer
    a += numpy.einsum('i,j->ij', x, y)
    return a

    def __init__(self, ndarray):

        self.metadata = _ndarray_meta(ndarray)
        self.size = ndarray.size * ndarray.itemsize
        # The following is guesswork
        # non contiguous fortran array (if ever such a thing exists)
        if numpy.isfortran(ndarray) and not ndarray.flags['F_CONTIGUOUS']:
            self.ndarray = numpy.asfortranarray(ndarray)
        # non contiguous C array
        elif not numpy.isfortran(ndarray) and not ndarray.flags['C_CONTIGUOUS']:
            self.ndarray = numpy.ascontiguousarray(ndarray)
        # contiguous fortran or C array, do nothing
        else:
            self.ndarray = ndarray
        self.ptr = self.ndarray.__array_interface__['data'][0]

 def verify_is_c_contiguous_and_is_not_fortran(x):
     if isinstance(x,CpuGpuArray):
         raise TypeError(type(x))            
     if np.isfortran(x):
         raise ValueError("Must be 'C' order")  
     if not x.flags.c_contiguous:
         raise ValueError("Must be 'C'-contiguous")      

    def tolocal(self, comm=MPI.COMM_WORLD):
        """Scatter a global image into local ones."""
        if hasattr(self, 'shape_global') and self.shape != self.shape_global or\
           hasattr(self, 'comm') and self.comm.Get_size() > 1:
            raise ValueError('This array is not a global image.')

        order = 'f' if np.isfortran(self) else 'c'
        if hasattr(self, 'empty'):
            output = self.empty(self.shape, dtype=self.dtype, order=order,
                                comm=comm)
        else:
            shape = split_shape(self.shape, comm)
            output = np.empty(shape, dtype=self.dtype, order=order)
            output = DistributedArray.__new__(DistributedArray, output,
                                              self.shape, comm)
        if hasattr(self, '__dict__'):
            for k,v in self.__dict__.items():
                setattr(output, k, v)
        output.comm = comm

        s = split_work(self.shape[0], comm=comm)
        n = s.stop - s.start
        output[0:n] = self[s.start:s.stop]
        if n < output.shape[0]:
            output[n:] = 0
        return output