def check_format(self, full_check=True):
        """check whether the matrix format is valid

        Parameters
        ==========

            - full_check : {bool}
                - True  - rigorous check, O(N) operations : default
                - False - basic check, O(1) operations

        """
        # use _swap to determine proper bounds
        major_name,minor_name = self._swap(('row','column'))
        major_dim,minor_dim = self._swap(self.shape)

        # index arrays should have integer data types
        if self.indptr.dtype.kind != 'i':
            warn("indptr array has non-integer dtype (%s)"
                    % self.indptr.dtype.name)
        if self.indices.dtype.kind != 'i':
            warn("indices array has non-integer dtype (%s)"
                    % self.indices.dtype.name)

        # only support 32-bit ints for now
        self.indptr = np.asarray(self.indptr, dtype=np.intc)
        self.indices = np.asarray(self.indices, dtype=np.intc)
        self.data = to_native(self.data)

        # check array shapes
        if np.rank(self.data) != 1 or np.rank(self.indices) != 1 or np.rank(self.indptr) != 1:
            raise ValueError('data, indices, and indptr should be rank 1')

        # check index pointer
        if (len(self.indptr) != major_dim + 1):
            raise ValueError("index pointer size (%d) should be (%d)" %
                                (len(self.indptr), major_dim + 1))
        if (self.indptr[0] != 0):
            raise ValueError("index pointer should start with 0")

        # check index and data arrays
        if (len(self.indices) != len(self.data)):
            raise ValueError("indices and data should have the same size")
        if (self.indptr[-1] > len(self.indices)):
            raise ValueError("Last value of index pointer should be less than "
                                "the size of index and data arrays")

        self.prune()

        if full_check:
            # check format validity (more expensive)
            if self.nnz > 0:
                if self.indices.max() >= minor_dim:
                    raise ValueError("%s index values must be < %d" %
                                        (minor_name,minor_dim))
                if self.indices.min() < 0:
                    raise ValueError("%s index values must be >= 0" %
                                        minor_name)
                if np.diff(self.indptr).min() < 0:
                    raise ValueError("index pointer values must form a "
                                        "non-decreasing sequence")

	def __init__(self,init_pos,init_measurement=[[]],init_weight=1,movement_weight=1,measurement_weight=1):
		if rank(init_pos)==0:
			init_pos=[init_pos]
		self.num_dim=len(init_pos)
		while (rank(init_measurement)<2):
			init_measurement=[init_measurement]
		if init_measurement==[[]]:
			self.num_landmarks=0
		else:
			self.num_landmarks=len(init_measurement[0])
		self.matrix_dim=self.num_dim+self.num_landmarks
		self.omega=zeros([self.matrix_dim,self.matrix_dim])
		self.xi=zeros([self.matrix_dim,1])
		self.init_weight=init_weight
		self.movement_weight=movement_weight
		self.measurement_weight=measurement_weight
		# initial position
		for i in range(self.num_dim):
			self.omega[i][i]+=self.init_weight
			self.xi[i][0]=init_pos[i]
			for j in range(self.num_landmarks):
				if init_measurement[i][j]!=None:
					self.omega[i][i]+=measurement_weight
					self.omega[i][self.num_dim+j]-=measurement_weight
					self.omega[self.num_dim+j][i]-=measurement_weight
					self.omega[self.num_dim+j][self.num_dim+j]+=measurement_weight
					self.xi[i][0]-=measurement_weight*init_measurement[i][j]
					self.xi[self.num_dim+j][0]+=measurement_weight*init_measurement[i][j]

def calc( absphot_vals_target, absphot_errs_target, absphot_vals_comparisons, absphot_errs_comparisons ):
    """
    Calculates relative fluxes given target and comparison flux time values. Also calculates the
    propagated formal uncertainty using the formal errors on each of the flux values.

    The required inputs are:
      ** absphot_vals_target - N-length np.array containing the absolute fluxes of the target star.
      ** absphot_errs_target - N-length np.array containing the error bars on the absolute fluxes of
         the target star.
      ** absphot_vals_comparisons - NxM np.array containing the absolute fluxes of M comparison stars.
      ** absphot_errs_comparisons - NxM np.array containing the error bars on the absolute fluxes of
         M comparison stars.

    All error calculation is done using standard uncertainty propagation assuming independence
    between points and quadrature sums.

    The function is defined in this way so that it's possible to have more flexibility in
    experimenting with different combinations of comparison stars elsewhere (eg. in other routines)
    if necessary.
    """
    
    # If there's only one comparison, just make sure everything is in the correct format:
    if ( np.rank( absphot_vals_comparisons )==1 ):
        absphot_vals_comparisons = np.reshape( absphot_vals_comparisons, [ len( absphot_vals_comparisons ), 1 ] )
    if ( np.rank( absphot_errs_comparisons )==1 ):
        absphot_errs_comparisons = np.reshape( absphot_errs_comparisons, [ len( absphot_errs_comparisons ), 1 ] )
    comparisons_sum_vals = np.sum( absphot_vals_comparisons, axis=1 )
    comparisons_sum_errs_sq = np.sum( absphot_errs_comparisons**2., axis=1 )
    comparisons_sum_errs = np.sqrt( comparisons_sum_errs_sq )
    relphot_vals = absphot_vals_target.flatten() / comparisons_sum_vals.flatten()
    relphot_errs = relphot_vals * np.sqrt( ( absphot_errs_target.flatten() / absphot_vals_target.flatten() )**2. \
                                        + ( comparisons_sum_errs.flatten() / comparisons_sum_vals.flatten() )**2. )
    return relphot_vals, relphot_errs

	def __init__(self,init_pos,init_meas,\
	init_weight=1,move_weight=1,meas_weight=1):
		"""
		sets up properties for the algorithm
		sets initial position and takes measurements there
		"""
		# allows for a variety of formatted input
		while rank(init_pos)<1:
			init_pos=[init_pos]
		while rank(init_meas)<1:
			init_meas=[init_meas]
		if rank(init_meas)==1:
			for i in range(len(init_meas)):
				if init_meas[i]!=None:
					init_meas[i]=[init_meas[i]]
		# set up properties for the algorithm
		self.num_dim=len(init_pos)
		self.num_landmarks=len(init_meas)
		self.matrix_dim=self.num_dim+self.num_landmarks
		self.omega=zeros([self.matrix_dim,self.matrix_dim])
		self.xi=zeros([self.matrix_dim,1])
		self.init_weight=init_weight
		self.move_weight=move_weight
		self.meas_weight=meas_weight
		# sets initial position
		self.setPosition(init_pos)
		# takes measurements at initial position
		self.measure(init_meas)

def orientation_product(T,Bb):
    """Computes the product of a tensor and a vector.

    Assumptions:
    None

    Source:
    N/A

    Inputs:
    T         [-] 3-dimensional array with rotation matrix
                  patterned along dimension zero
    Bb        [-] 3-dimensional vector

    Outputs:
    C         [-] transformed vector

    Properties Used:
    N/A
    """            
    
    assert np.rank(T) == 3
    
    if np.rank(Bb) == 3:
        C = np.einsum('aij,ajk->aik', T, Bb )
    elif np.rank(Bb) == 2:
        C = np.einsum('aij,aj->ai', T, Bb )
    else:
        raise Exception , 'bad B rank'
        
    return C

def coordinate_transform(old_coord_data, iapp=None):
	# Translate, then transform
	assert (type(iapp) is float or type(iapp) is np.float64 or iapp == None)
	saddle = saddle_point(iapp)
	eigenvectors = get_eigenvectors(saddle, iapp)

	eigenvectors[:,0] = eigenvectors[:,0]/linalg.norm(eigenvectors[:,0])
	eigenvectors[:,1] = eigenvectors[:,1]/linalg.norm(eigenvectors[:,1])

	raw_matrix = eigenvectors
	transformed = linalg.inv(raw_matrix)

	if len(old_coord_data)==2 and np.rank(old_coord_data)==1:
		new_coord = np.zeros([2,1])
		new_coord[0] = old_coord_data[0] - saddle[0]
		new_coord[1] = old_coord_data[1] - saddle[1]
		new_coord_final = np.dot(transformed, np.array([new_coord[0], new_coord[1]]))
		return new_coord_final


	elif old_coord_data.shape[0] == 2 and np.rank(old_coord_data)==2:
		old_coord_data[0] = old_coord_data[0] - saddle[0]
		old_coord_data[1] = old_coord_data[1] - saddle[1]
		for i in range(len(old_coord_data[0])):
			old_coord_data[:,i] = np.dot(transformed, old_coord_data[:,i])
		return old_coord_data

	elif old_coord_data.shape[1] == 2 and np.rank(old_coord_data)==2:
		old_coord_data[:,0] = old_coord_data[:,0] - saddle[0]
		old_coord_data[:,1] = old_coord_data[:,1] - saddle[1]
		for i in range(len(old_coord_data)):
			old_coord_data[i] = np.dot(transformed, old_coord_data[i])
		return old_coord_data
	else:
		raise Exception("couldn't find appropriate manipulation for", old_coord_data.shape)

def find_duplicates_others(source_array, other_array):
    """find indices of duplicate values in src_array against other_array
    source_array - numpy array to be checked
    other_array - numpy array to be checked againt, must have compatiable shape[0]
                   with src_array
    """

    if other_array is None or len(other_array) == 0:
        return zeros(source_array.shape, dtype="int32")

    if other_array.shape[0] <> source_array.shape[0]:
        raise ValueError, "Arrays have incompatible shapes"
    source_array_rank = rank(source_array)
    if rank(source_array) < rank(other_array):
        src_array = source_array[:, newaxis]
        oth_array = other_array
    elif rank(source_array) > rank(other_array):
        src_array = source_array
        oth_array = other_array[:, newaxis]

    is_duplicates = equal(src_array, oth_array)
    duplicate_indicator = sum(is_duplicates, axis=1)

    #    if duplicate_indicator.ndim > source_array_rank:
    #        reshape(duplicate_indicator, shape=(src_array.size,))
    return duplicate_indicator

def get_vel(pos, time, tax=0, spax=-1):
    ''' Get instantaneous velocity

    Parameters
    ----------
    time : array_like

    pos : array_like

    tax : int, optional
      time axis, defaults to 0
      has to be suitable for both pos and time,
      so probably needs to be +ve, i.e. indexed from beginning
    spax : int, optional
      space axis in pos, defaults to -1, i.e. last axis
    '''
    dp = np.diff(pos, axis=tax)
    dt = np.diff(time, axis=tax)

    if np.rank(dp) != np.rank(dt):
        if spax < 0:
            spax = len(pos.shape) + spax
        dts = [slice(None) for x in dt.shape]
        dts.insert(spax, None)
    else:
        dts = [slice(None) for x in dt.shape]
    return dp / dt[dts]

def _format_for_gam2(count, time, pos):
    assert np.rank(count) == np.rank(time) == (np.rank(pos) - 1)
    assert count.shape[0] == time.shape[0] == pos.shape[0]
    assert (count.shape[1] + 1) == time.shape[1] == pos.shape[1]

    y = count.flatten()[:,None]
    npt = y.size
    tax, spax = 1, -1
    
    # don't use real time as won't compare equivalent portions of trials
    # t  = edge2cen(time, axis=tax) # (ntrial, nbin)
    # subtract offset to get all relative times
    #t = (t - t[:,:1]).flatten()[:,None]
    
    # instead use bin numbers
    ntrial, nbin = count.shape
    t = np.tile(np.arange(nbin, dtype=float)[None], (ntrial, 1))
    t = t.flatten()[:,None]

    # also create trial numbers for pulling out appropriate later on
    # these don't correspond to original trial numbers because original trials
    # have been permuted before getting here
    tr = np.tile(np.arange(ntrial), (nbin, 1)).T.flatten()
    
    d  = kin.get_dir(pos, tax=tax, spax=spax).reshape(npt, 3)
    p  = edge2cen(pos, axis=tax).reshape(npt, 3)
    v  = kin.get_vel(pos, time, tax=tax, spax=spax).reshape(npt, 3)
    sp = kin.get_speed(pos, time, tax=tax, spax=spax).flatten()[:,None]

    # q is a second direction set-of-columns for deviance calculation
    q  = kin.get_dir(pos, tax=tax, spax=spax).reshape(npt, 3)
    return np.concatenate([y,t,d,p,v,sp,q], axis=1), tr

def coordinate_reverse_transform(old_coord_data, iapp):
	# Translate, then transform
	transformed = np.array([[1.245218544034379,-0.7978622303790213],[-0.3779055518316757,1.4298048589659018]])
	transformed = linalg.inv(transformed)

	if len(old_coord_data)==2 and np.rank(old_coord_data)==1:
		old_coord_data = np.dot(transformed, old_coord_data)
		old_coord_data[0] = old_coord_data[0] + 0.4 - 2*iapp
		old_coord_data[1] = old_coord_data[1] + 0.15 - iapp
		return old_coord_data

	elif old_coord_data.shape[1] == 2 and np.rank(old_coord_data)==2:
		for i in range(len(old_coord_data)):
			old_coord_data[i] = np.dot(transformed, old_coord_data[i])
		old_coord_data[:,0] = old_coord_data[:,0] + 0.4 - 2*iapp
		old_coord_data[:,1] = old_coord_data[:,1] + 0.15 - iapp
		return old_coord_data

	elif old_coord_data.shape[0] == 2 and np.rank(old_coord_data)==2:
		for i in range(len(old_coord_data[0])):
			old_coord_data[:,i] = np.dot(transformed, old_coord_data[:,i])
		old_coord_data[0] = old_coord_data[0] + 0.4 - 2*iapp
		old_coord_data[1] = old_coord_data[1] + 0.15 - iapp
		return old_coord_data

	else:
		raise Exception("couldn't find appropriate manipulation for", old_coord_data.shape)

def coordinate_transform(old_coord_data, iapp):
	# Translate, then transform
	transformed = np.array([[1.245218544034379,-0.7978622303790213],[-0.3779055518316757,1.4298048589659018]])

	if len(old_coord_data)==2 and np.rank(old_coord_data)==1:
		new_coord = np.zeros([2,1])
		#print new_coord
		#print old_coord_data
		new_coord[0] = old_coord_data[0] - 0.4 + 2*iapp
		new_coord[1] = old_coord_data[1] - 0.15 + iapp
		new_coord_final = np.dot(transformed, np.array([new_coord[0], new_coord[1]]))
		return new_coord_final

	elif old_coord_data.shape[1] == 2 and np.rank(old_coord_data)==2:
		old_coord_data[:,0] = old_coord_data[:,0] - 0.4 + 2*iapp
		old_coord_data[:,1] = old_coord_data[:,1] - 0.15 + iapp
		for i in range(len(old_coord_data)):
			old_coord_data[i] = np.dot(transformed, old_coord_data[i])
		return old_coord_data

	elif old_coord_data.shape[0] == 2 and np.rank(old_coord_data)==2:
		old_coord_data[0] = old_coord_data[0] - 0.4 + 2*iapp
		old_coord_data[1] = old_coord_data[1] - 0.15 + iapp
		for i in range(len(old_coord_data[0])):
			old_coord_data[:,i] = np.dot(transformed, old_coord_data[:,i])
		return old_coord_data			

	else:
		raise Exception("couldn't find appropriate manipulation for", old_coord_data.shape)

    def getnnz(self, axis=None):
        """Get the count of explicitly-stored values (nonzeros)

        Parameters
        ----------
        axis : None, 0, or 1
            Select between the number of values across the whole matrix, in
            each column, or in each row.
        """
        if axis is None:
            nnz = len(self.data)
            if nnz != len(self.row) or nnz != len(self.col):
                raise ValueError('row, column, and data array must all be the '
                                 'same length')

            if np.rank(self.data) != 1 or np.rank(self.row) != 1 or \
               np.rank(self.col) != 1:
                raise ValueError('row, column, and data arrays must have '
                                 'rank 1')

            return int(nnz)

        if axis < 0:
            axis += 2
        if axis == 0:
            return _compat_bincount(downcast_intp_index(self.col),
                                    minlength=self.shape[1])
        elif axis == 1:
            return _compat_bincount(downcast_intp_index(self.row),
                                    minlength=self.shape[0])
        else:
            raise ValueError('axis out of bounds')

def project_scores_to_var_space(score, weight, data):
    '''
    Project reduced scores, via reduced weights, up to neuron space    
    
    Parameters
    ----------
    score : ndarray
      shape (npc, nobs), i.e. (nscore, ntask [* nrep] * nbin)
    weight : ndarray
      shape (npc, nvar), i.e. (nscore, nunit)
    data : ndarray
      shape (nvar, nobs), data from which to get mean
      
    Returns
    -------
    projected : ndarray
      shape (nvar, nobs), i.e. (nunit, ntask * nbin)
    '''
    assert np.rank(score) == np.rank(weight) == np.rank(data) == 2
    assert score.shape[0] == weight.shape[0] # npc
    assert score.shape[1] == data.shape[1]   # nobs
    assert weight.shape[1] == data.shape[0]  # nvar
    
    # take average over observations
    mean = stats.nanmean(data, axis=1)
    return (np.dot(weight.T, score) + mean[:,None])

def cube_array_search(k_face_array,k_faces):
    """
    Find the row indices (of s) corresponding to the
    cubes stored in the rows of cube array v.
    It is assumed that the rows of s are sorted in
    lexicographical order.

    Example:

      k_face_array = array([[0,0,0],[0,0,1],[0,1,0],[1,0,1]])
      k_faces = array([[0,1,0],[0,0,1]])
      cube_array_searchsorted(k_face_array,k_faces)

    Returns:

      array([2,1])

    """
    if rank(k_face_array) != 2 or rank(k_faces) != 2:
        raise ValueError,'expected rank 2 arrays'

    if k_face_array.shape[1] != k_faces.shape[1]:
        raise ValueError,'number of columns must agree'

    # a dense array used to lookup k_face_array row indices 
    lookup_grid_dimensions = k_face_array.max(axis=0) + 1
    
    lookup_grid = empty(lookup_grid_dimensions,dtype=k_faces.dtype)
    lookup_grid[:] = -1
    lookup_grid[hsplit(k_face_array,k_face_array.shape[1])] = arange(k_face_array.shape[0],dtype=k_faces.dtype).reshape((-1,1))
    row_indices = lookup_grid[hsplit(k_faces,k_faces.shape[1])].reshape((-1))

    return row_indices

    def check_format(self, full_check=True):
        """check whether the matrix format is valid

            *Parameters*:
                full_check:
                    True  - rigorous check, O(N) operations : default
                    False - basic check, O(1) operations

        """
        M,N = self.shape
        R,C = self.blocksize

        # index arrays should have integer data types
        if self.indptr.dtype.kind != 'i':
            warn("indptr array has non-integer dtype (%s)" \
                    % self.indptr.dtype.name )
        if self.indices.dtype.kind != 'i':
            warn("indices array has non-integer dtype (%s)" \
                    % self.indices.dtype.name )

        # only support 32-bit ints for now
        self.indptr  = np.asarray(self.indptr, np.intc)
        self.indices = np.asarray(self.indices, np.intc)
        self.data    = to_native(self.data)

        # check array shapes
        if np.rank(self.indices) != 1 or np.rank(self.indptr) != 1:
            raise ValueError,"indices, and indptr should be rank 1"
        if np.rank(self.data) != 3:
            raise ValueError,"data should be rank 3"

        # check index pointer
        if (len(self.indptr) != M/R + 1 ):
            raise ValueError, \
                "index pointer size (%d) should be (%d)" % \
                 (len(self.indptr), M/R + 1)
        if (self.indptr[0] != 0):
            raise ValueError,"index pointer should start with 0"

        # check index and data arrays
        if (len(self.indices) != len(self.data)):
            raise ValueError,"indices and data should have the same size"
        if (self.indptr[-1] > len(self.indices)):
            raise ValueError, \
                  "Last value of index pointer should be less than "\
                  "the size of index and data arrays"

        self.prune()

        if full_check:
            #check format validity (more expensive)
            if self.nnz > 0:
                if self.indices.max() >= N/C:
                    print "max index",self.indices.max()
                    raise ValueError, "column index values must be < %d" % (N/C)
                if self.indices.min() < 0:
                    raise ValueError, "column index values must be >= 0"
                if diff(self.indptr).min() < 0:
                    raise ValueError,'index pointer values must form a " \

    def SI(self):
        """
        Return the quantity in SI unit. If the quantity has
        no units, the quantity itself is returned, otherwise, a
        shallow copy is returned.

        Example
        -------
        >>> print(Quantity(1., 'km').SI)
        1000.0 m
        """
        if len(self._unit) == 0:
            return self

        factor = Quantity(1., '')
        for key, val in self._unit.items():

            # check if the unit is a local derived unit
            newfactor = _check_du(self, key, val, self.derived_units)

            # check if the unit is a global derived unit
            if newfactor is None:
                newfactor = _check_du(self, key, val, units)
                
            # if the unit is not derived, we add it to the dictionary
            if newfactor is None:            
                _multiply_unit_inplace(factor._unit, key, val)
                continue

            # factor may be broadcast
            factor = factor * newfactor

        if np.rank(self) == 0 or np.rank(factor) == 0:
            result = self * factor.magnitude
            result._unit = factor._unit
            return result

        # make sure that the unit factor can be broadcast
        if any([d1 not in (1,d2)
                for (d1,d2) in zip(self.shape[0:factor.ndim], factor.shape)]):
            raise ValueError("The derived unit '" + key + "' has a shape '" + \
                str(newfactor.shape) + "' which is incompatible with the dime" \
                "nsion(s) along the first axes '" + str(self.shape) + "'.")

        # Python's broadcast operates by adding slow dimensions. Since we
        # need to use a unique conversion factor along the fast dimensions
        # (ex: second axis in an array tod[detector,time]), we need to perform
        # our own broadcast by adding fast dimensions
        if np.any(self.shape[0:factor.ndim] != factor.shape):
            broadcast = np.hstack([factor.shape,
                                  np.ones(self.ndim-factor.ndim,int)])
            result = self * np.ones(broadcast)
        else:
            result = self.copy()
        result.T[:] *= factor.magnitude.T
        result._unit = factor._unit

        return result

def fftconvolve(
    in1,
    in2,
    mode="full",
    preprocessed_input=None
    ):
    """Convolve two N-dimensional arrays using FFT.

    Very similar to scipy.signal.fftconvolve, but with some performance 
    customizations to the padding, and optional preprocessing.

    scipy.signal.fftconvolve currently pads to a power of 2 for fast
    FFTs; this is often suboptimal, since numpy's FFTs are fast for
    array sizes that are products of small primes. Our version tries
    to guess a good padding size.

    Often, we find ourselves convolving many different in1's with the
    same in2. The 'preprocessed_input' argument lets us save
    computation time by passing a preprocessed version of in2.
    """
    in1 = np.asarray(in1)
    in2 = np.asarray(in2)

    if np.rank(in1) == np.rank(in2) == 0:  # scalar inputs
        return in1 * in2
    elif not in1.ndim == in2.ndim:
        raise ValueError("in1 and in2 should have the same rank")
    elif in1.size == 0 or in2.size == 0:  # empty arrays
        return array([])

    if preprocessed_input is None:
        preprocessed_input = preprocess_input_2(in1, in2)
    (in2_preprocessed,
     in2_preprocessed_is_complex,
     size,
     fsize,
     s1,
     s2
     ) = preprocessed_input

    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      in2_preprocessed_is_complex)

    fslice = tuple([slice(0, int(sz)) for sz in size])
    if not complex_result:
        ret = np.fft.irfftn(np.fft.rfftn(in1, fsize) * in2_preprocessed, fsize
                            )[fslice].copy()
        ret = ret.real
    else:
        ret = np.fft.ifftn(np.fft.fftn(in1, fsize) * in2_preprocessed
                           )[fslice].copy()

    if mode == "full":
        return ret
    elif mode == "same":
        return _centered(ret, s1)
    elif mode == "valid":
        return _centered(ret, s1 - s2 + 1)

def add_column(cur_vars, more_vars):
    if np.rank(more_vars) == 1:
        more_vars.shape = list(more_vars.shape) + [1,]
    if cur_vars != None: # must have same number of tasks
        assert np.rank(more_vars) == 2
        assert cur_vars.shape[0] == more_vars.shape[0]
        return np.concatenate((cur_vars, more_vars), axis=-1)
    else:
        return more_vars

    def getnnz(self):
        nnz = len(self.data)
        if nnz != len(self.row) or nnz != len(self.col):
            raise ValueError('row, column, and data array must all be the same length')

        if np.rank(self.data) != 1 or np.rank(self.row) != 1 or np.rank(self.col) != 1:
            raise ValueError('row, column, and data arrays must have rank 1')

        return int(nnz)

    def check_format(self, full_check=True):
        """check whether the matrix format is valid

            *Parameters*:
                full_check:
                    True  - rigorous check, O(N) operations : default
                    False - basic check, O(1) operations

        """
        M,N = self.shape
        R,C = self.blocksize

        # index arrays should have integer data types
        if self.indptr.dtype.kind != 'i':
            warn("indptr array has non-integer dtype (%s)"
                    % self.indptr.dtype.name)
        if self.indices.dtype.kind != 'i':
            warn("indices array has non-integer dtype (%s)"
                    % self.indices.dtype.name)

        idx_dtype = get_index_dtype((self.indices, self.indptr))
        self.indptr = np.asarray(self.indptr, dtype=idx_dtype)
        self.indices = np.asarray(self.indices, dtype=idx_dtype)
        self.data = to_native(self.data)

        # check array shapes
        if np.rank(self.indices) != 1 or np.rank(self.indptr) != 1:
            raise ValueError("indices, and indptr should be rank 1")
        if np.rank(self.data) != 3:
            raise ValueError("data should be rank 3")

        # check index pointer
        if (len(self.indptr) != M//R + 1):
            raise ValueError("index pointer size (%d) should be (%d)" %
                                (len(self.indptr), M//R + 1))
        if (self.indptr[0] != 0):
            raise ValueError("index pointer should start with 0")

        # check index and data arrays
        if (len(self.indices) != len(self.data)):
            raise ValueError("indices and data should have the same size")
        if (self.indptr[-1] > len(self.indices)):
            raise ValueError("Last value of index pointer should be less than "
                                "the size of index and data arrays")

        self.prune()

        if full_check:
            # check format validity (more expensive)
            if self.nnz > 0:
                if self.indices.max() >= N//C:
                    raise ValueError("column index values must be < %d (now max %d)" % (N//C, self.indices.max()))
                if self.indices.min() < 0:
                    raise ValueError("column index values must be >= 0")
                if np.diff(self.indptr).min() < 0:
                    raise ValueError("index pointer values must form a "
                                        "non-decreasing sequence")