def prob_of_label(vol, labelvol):
    """
    compute the probability of the labels in labelvol in each of the volumes in vols

    Parameters:
        vol (float numpy array of dim (nd + 1): volume with a prob dist at each voxel in a nd vols
        labelvol (int numpy array of dim nd): nd volume of labels

    Returns:
        nd volume of probabilities
    """

    # check dimensions
    ndims = np.ndim(labelvol)
    assert np.ndim(vol) == ndims + 1, "vol dimensions do not match [%d] vs [%d]" % (np.ndim(vol)-1, ndims)
    shp = vol.shape
    nb_voxels = np.prod(shp[0:ndims])
    nb_labels = shp[-1]

    # reshape volume to be [nb_voxels, nb_labels]
    flat_vol = np.reshape(vol, (nb_voxels, nb_labels))

    # normalize accross second dimension
    rows_sums = flat_vol.sum(axis=1)
    flat_vol_norm = flat_vol / rows_sums[:, np.newaxis]

    # index into the flattened volume
    idx = list(range(nb_voxels))
    v = flat_vol_norm[idx, labelvol.flat]
    return np.reshape(v, labelvol.shape)

def py_vq(obs, code_book):
    """ Python version of vq algorithm.

    The algorithm computes the euclidian distance between each
    observation and every frame in the code_book.

    Parameters
    ----------
    obs : ndarray
        Expects a rank 2 array. Each row is one observation.
    code_book : ndarray
        Code book to use. Same format than obs. Should have same number of
        features (eg columns) than obs.

    Returns
    -------
    code : ndarray
        code[i] gives the label of the ith obversation, that its code is
        code_book[code[i]].
    mind_dist : ndarray
        min_dist[i] gives the distance between the ith observation and its
        corresponding code.

    Notes
    -----
    This function is slower than the C version but works for
    all input types.  If the inputs have the wrong types for the
    C versions of the function, this one is called as a last resort.

    It is about 20 times slower than the C version.

    """
    # n = number of observations
    # d = number of features
    if np.ndim(obs) == 1:
        if not np.ndim(obs) == np.ndim(code_book):
            raise ValueError(
                    "Observation and code_book should have the same rank")
        else:
            return _py_vq_1d(obs, code_book)
    else:
        (n, d) = shape(obs)

    # code books and observations should have same number of features and same
    # shape
    if not np.ndim(obs) == np.ndim(code_book):
        raise ValueError("Observation and code_book should have the same rank")
    elif not d == code_book.shape[1]:
        raise ValueError("Code book(%d) and obs(%d) should have the same " \
                         "number of features (eg columns)""" %
                         (code_book.shape[1], d))

    code = zeros(n, dtype=int)
    min_dist = zeros(n)
    for i in range(n):
        dist = np.sum((obs[i] - code_book) ** 2, 1)
        code[i] = argmin(dist)
        min_dist[i] = dist[code[i]]

    return code, sqrt(min_dist)

def _test():
    imdata = np.array([[1,2,3],[4,5,6]])
    immask = imdata>2
    print imdata
    print immask
    print imdata[immask]
    print np.ndim(immask)

def GetLaneFromStdDeviation(cluster):

    if cluster is None:
        print "GetLaneFromMedian: cluster is empty"
        return None

    cluster_removed_outliers = RemoveOutliers(cluster)

    if cluster_removed_outliers is None:
        # IN THIS CASE NO VALUES ARE WITHIN ONE STD DEVIATION!
        if np.ndim(cluster) is 1:
            distances = cluster[0]
            angles = cluster[1]
        else:
            distances = cluster[:, 0]
            angles = cluster[:, 1]
        # return None
        # print 'Aww'
    elif np.ndim(cluster_removed_outliers) is 1:
        distances = cluster_removed_outliers[0]
        angles = cluster_removed_outliers[1]
    else:
        distances = cluster_removed_outliers[:, 0]
        angles = cluster_removed_outliers[:, 1]

    avg_angle = np.mean(angles)
    avg_distance = np.mean(distances)

    return np.array([avg_distance, avg_angle])

def as2d(a):
  if np.ndim(a) == 0:
    return a[np.newaxis,np.newaxis]
  elif np.ndim(a) == 1:
    return a[:,np.newaxis]
  else:
    return a

def _populate_vars(dz, zdr, kdp, rho, ldr, T, verbose):

    """
    Check for presence of each var, and update dicts as needed.
    Flattens multi-dimensional arrays to optimize processing.
    The output array from csu_fhc_summer() will be re-dimensionalized later.
    """
    varlist = [dz, zdr, kdp, rho, ldr, T]
    keylist = ['DZ', 'DR', 'KD', 'RH', 'LD', 'T']
    fhc_vars = {}
    radar_data = {}
    for i, key in enumerate(keylist):
        var = varlist[i]
        if var is not None:
            if key == 'DZ':
                shp = np.shape(var)
                sz = np.size(var)
            if np.ndim(var) > 1:
                radar_data[key] = np.array(var).ravel().astype('float32')
            elif np.ndim(var) == 1:
                radar_data[key] = np.array(var).astype('float32')
            else:
                radar_data[key] = np.array([var]).astype('float32')
            fhc_vars[key] = 1
        else:
            fhc_vars[key] = 0
    if verbose:
        print('USING VARIABLES: ', fhc_vars)
    return radar_data, fhc_vars, shp, sz

def nan_detrend(x, y, deg=1):
    """Subtract a polynomial fit from the data, ignoring NaNs.

    Parameters
    ----------
    x : array_like
        x data.
    y : array_like
        Data to detrend.
    deg : int
        Degree of polynomial to subtract. (Can be zero i.e. constant)


    Returns
    -------
    y_out : numpy.array
        Detrended data.

    """
    y_out = np.nan*np.zeros_like(y)

    if np.ndim(x) == 1:
        nans = np.isnan(x) | np.isnan(y)
        p = nan_polyfit(x, y, deg)
        y_out[~nans] = y[~nans] - np.polyval(p, x[~nans])
    elif np.ndim(x) == 2:
        for i in xrange(x.shape[1]):
            nans = np.isnan(x[:, i]) | np.isnan(y[:, i])
            p = nan_polyfit(x[:, i], y[:, i], deg)
            y_out[~nans, i] = y[~nans, i] - np.polyval(p, x[~nans, i])
    else:
        raise RuntimeError('Arguments must be 1 or 2 dimensional arrays.')

    return y_out

def read_h5_stack(fn, group='stack', crop=[None]*6, **kwargs):
    """Read a volume in HDF5 format into numpy.ndarray.

    Parameters
    ----------
    fn : string
        The filename of the input HDF5 file.
    group : string, optional (default 'stack')
        The group within the HDF5 file containing the dataset.
    crop : list of int, optional (default '[None]*6', no crop)
        A crop to get of the volume of interest. Only available for 2D and 3D
        volumes.

    Returns
    -------
    stack : numpy ndarray
        The stack contained in fn, possibly cropped.
    """
    fn = os.path.expanduser(fn)
    dset = h5py.File(fn, 'r')
    if group not in dset:
        raise ValueError("HDF5 file (%s) doesn't have group (%s)!" % 
                            (fn, group))
    a = dset[group]
    if ndim(a) == 2:
        xmin, xmax, ymin, ymax = crop[:4]
        a = a[xmin:xmax, ymin:ymax]
    elif ndim(a) == 3:
        xmin, xmax, ymin, ymax, zmin, zmax = crop
        a = a[xmin:xmax, ymin:ymax, zmin:zmax]
    stack = array(a)
    dset.close()
    return stack

    def project(self, *args, **kwargs):
        """ Project molecule.

        Parameters
        ----------
        mol : :class:`Molecule <htmd.molecule.molecule.Molecule>`
            A :class:`Molecule <htmd.molecule.molecule.Molecule>` object to project.

        Returns
        -------
        data : np.ndarray
            An array containing the projected data.
        """
        # -------------- DEPRECATION PATCH --------------
        if isinstance(self, np.ndarray) or isinstance(self, Molecule):
            from warnings import warn
            warn('Static use of the project method will be deprecated in the next version of HTMD. '
                 'Please change your projection code according to the tutorials on www.htmd.org')
            data = _MetricDistanceOld.project(self, *args, **kwargs)
            logger.warning('Static use of the project method will be deprecated in the next version of HTMD. '
                            'Please change your projection code according to the tutorials on www.htmd.org')
            return data
        # ------------------ CUT HERE -------------------
        mol = args[0]
        (sel1, sel2) = self._getSelections(mol)

        if np.ndim(sel1) == 1 and np.ndim(sel2) == 1:  # normal distances
            metric = pp_calcDistances(mol, sel1, sel2, self.metric, self.threshold, self.pbc, truncate=self.truncate)
        else:  # minimum distances by groups
            metric = pp_calcMinDistances(mol, sel1, sel2)

        return metric

def sample_vMF(theta, kappa,size=1):
    """
      Sampling from vMF
      This is based on the implementation I found online here:
          http://stats.stackexchange.com/questions/156729/sampling-from-von-mises-fisher-distribution-in-python
          (**** NOTE THE FIX BY KEVIN *****)
      which is based on : 
            Directional Statistics (Mardia and Jupp, 1999) and on the Ulrich-Wood's algorithm for sampling. 


    """

    warn('Not sure about sampling vMF, use with caution!!!! ')
    #print "kappa : ", kappa
    #print "norm direction :" , np.linalg.norm(theta)
    np.testing.assert_array_almost_equal( np.linalg.norm(theta) , 1 )
    #print "kappa : ", kappa
    assert kappa > 0 , "kappa must be positive !"

    if np.ndim(theta)==2:
       theta = np.squeeze(theta)
    assert np.ndim(theta)==1, "theta should be one one-dimensional!"
    
    res_sampling = _rvMF(size, kappa * theta)

    return np.vstack(res_sampling)    

 def helper():
     flist = glob.iglob(os.path.join(directory, "*.img"))
     flist = list(itertools.ifilter(lambda f: os.path.getmtime(f) > time0[0], 
                                    flist))
     
     if not flist:
         return
     flist.sort(key=os.path.getmtime)
     
     fpath = flist[-1]
     time0[0] = os.path.getmtime(fpath)
     data = open_img(fpath).data
     data = np.rot90(data, -1)
     if np.ndim(iw.image)==2 and iw.image.shape==data.shape:
         newdata = np.concatenate([data[np.newaxis], iw.image[np.newaxis]], 
                                  axis=0)
     elif np.ndim(iw.image)==3 and iw.image.shape[1:3]==data.shape:
         newdata = np.concatenate([data[np.newaxis], iw.image], axis=0)
     elif np.ndim(data)==2:
         newdata = data
     else:
         return
     print fpath
     iw.setImage(newdata, autoLevels=False)
     iw.labels.insert(0, fpath)
     iw.win.setWindowTitle(fpath)

def fastdtw(x, y, dist):
    """
    Computes Dynamic Time Warping (DTW) of two sequences in a faster way.
    Instead of iterating through each element and calculating each distance,
    this uses the cdist function from scipy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html)
    :param array x: N1*M array
    :param array y: N2*M array
    :param string or func dist: distance parameter for cdist. When string is given, cdist uses optimized functions for the distance metrics.
    If a string is passed, the distance function can be 'braycurtis', 'canberra', 'chebyshev', 'cityblock', 'correlation', 'cosine', 'dice', 'euclidean', 'hamming', 'jaccard', 'kulsinski', 'mahalanobis', 'matching', 'minkowski', 'rogerstanimoto', 'russellrao', 'seuclidean', 'sokalmichener', 'sokalsneath', 'sqeuclidean', 'wminkowski', 'yule'.
    Returns the minimum distance, the cost matrix, the accumulated cost matrix, and the wrap path.
    """
    assert len(x)
    assert len(y)
    if ndim(x)==1:
        x = x.reshape(-1,1)
    if ndim(y)==1:
        y = y.reshape(-1,1)
    r, c = len(x), len(y)
    D0 = zeros((r + 1, c + 1))
    D0[0, 1:] = inf
    D0[1:, 0] = inf
    D1 = D0[1:, 1:]
    D0[1:,1:] = cdist(x,y,dist)
    C = D1.copy()
    for i in range(r):
        for j in range(c):
            D1[i, j] += min(D0[i, j], D0[i, j+1], D0[i+1, j])
    if len(x)==1:
        path = zeros(len(y)), range(len(y))
    elif len(y) == 1:
        path = range(len(x)), zeros(len(x))
    else:
        path = _traceback(D0)
    return D1[-1, -1] / sum(D1.shape), C, D1, path

def sparse_jac_repeat(tape_tag, x, nnz, rind, cind, values):
    """
    computes sparse Jacobian J for a function F:R^N -> R^M with
    the sparsity pattern that has been computed previously (e.g. by calling sparse_jac_no_repeat)

    I guess it also reuses the options that have been set previously. So it would be not necessary to set the options again.

    
    [nnz, rind, cind, values] = sparse_jac_repeat(tape_tag, x, rind, cind, values)

    INPUT:
    The base point x at which the Jacobian should be computed, i.e. J = dF(x)/dx
    
    OUTPUT:
    nnz is the guessed number of nonzeros in the Jacobian. This can be larger than the true number of nonzeros.

    sparse matrix representation in standard format:
    rind is an nnz-array of row indices
    cind is an nnz-array of column indices
    values are the corresponding Jacobian entries
    """
    assert type(tape_tag) == int
    assert type(nnz) == int

    assert numpy.ndim(x) == 1
    assert numpy.ndim(rind) == 1
    assert numpy.ndim(cind) == 1
    assert numpy.ndim(values) == 1

    x = numpy.asarray(x, dtype=float)
    rind = numpy.asarray(rind, dtype=numpy.uint32)
    cind = numpy.asarray(cind, dtype=numpy.uint32)
    values = numpy.asarray(values, dtype=float)

    return _colpack.sparse_jac_repeat(tape_tag, x, nnz, rind, cind, values)

def interpolate_1x1(lon, lat, data, lon0=0.):
    """ interpolate on a standard grid
    """
    # already standard
    if is_standard_grid(lon, lat, lon0=lon0):
        return lon, lat, data

    # 
    lon, data = rectify_longitude_data(lon, data, lon0)

    res = 1
    nx = int(360/res)
    ny = int(180/res)

    lon360 = np.linspace(lon0+0.5,lon0+359.5,nx)
    lat180 = np.linspace(-89.5,89.5,ny)

    print "Interpolate onto a standard 1deg grid...",
    if np.ndim(data) == 2:
        if isinstance(data, np.ma.MaskedArray):
            data = data.filled(np.nan)
        data = interp2(lon, lat, data, lon360,lat180, order=0)
        lon, lat = lon360, lat180

    elif np.ndim(data) == 3:
        nt = np.size(data,0)
        data180x360 = np.zeros((nt,ny,nx))
        for i in range(nt):
            data180x360[i,:,:] = interp2(lon, lat, np.squeeze(data[i,:,:]), lon360,lat180)
        lon, lat, data = lon360, lat180, data180x360
    print "Done"

    return lon, lat, data

    def getMapping(self, mol):
        (sel1, sel2) = self._getSelections(mol)

        if np.ndim(sel1) == 2:
            protatoms = np.ones(len(sel1)) * -1
            for i in range(np.size(sel1, 1)):
                protatoms[i] = np.where(sel1[i] == True)[0][0]
        else:
            protatoms = np.where(sel1)[0]
        if np.ndim(sel2) == 2:
            ligatoms = np.ones(len(sel1)) * -1
            for i in range(np.size(sel2, 1)):
                ligatoms[i] = np.where(sel2[i] == True)[0][0]
        else:
            ligatoms = np.where(sel2)[0]

        numatoms1 = len(protatoms)
        numatoms2 = len(ligatoms)

        if np.array_equal(protatoms, ligatoms):
            map = np.zeros((numatoms1 * (numatoms1-1) / 2, 2), dtype=int)
            start = 0
            for i in range(numatoms1):
                finish = start + numatoms1 - i - 1
                map[start:finish, 0] = protatoms[i]
                map[start:finish, 1] = protatoms[i+1:]
                start = finish
        else:
            map = np.zeros((numatoms1 * numatoms2, 2), dtype=int)
            for i in range(numatoms2):
                start = i * numatoms1
                finish = (i+1) * numatoms1
                map[start:finish, 0] = protatoms
                map[start:finish, 1] = ligatoms[i]
        return map

def xcorr(x,y,**kwargs):
	"""cross correlation by rfft"""
	x = np.asarray(x)
	y = np.asarray(y)
	if np.ndim(x) == np.ndim(y):
		shape=kwargs.get('shape',np.max((x.shape, y.shape), axis = 0))
		return np.fft.irfftn(np.conjugate(np.fft.rfftn(x,s=shape))*np.fft.rfftn(y,s=shape))
	elif np.ndim(y) == 1:
		axis = kwargs.get('axis', 0)
		shape=kwargs.get('shape', max(x.shape[axis], len(y)))
		shape+=shape%2
		outshape = np.array(x.shape[:])
		outshape[axis] = shape
		out = np.zeros(outshape)
		y = np.fft.ifftshift(np.pad(y, pad_width = (int((shape-len(y)+1)/2), int((shape-len(y))/2)), mode = 'constant'))
		y_fft = np.fft.rfft(y, n=shape)
		x_fft = np.fft.rfft(x, n=shape, axis=axis)
		if axis == 0:
			for ii in range(len(x_fft[0])):
				out[:,ii] = np.fft.irfft(x_fft[:,ii]*np.conjugate(y_fft))
		else:
			for ii in range(len(x_fft)):
				out[ii] = np.fft.irfft(x_fft[ii]*np.conjugate(y_fft))
		return out
	else:
		raise ValueError('Only inputs with dimensions of 1 or 2 can be processed.')

    def transformImage(self, containerImage):
            if type(containerImage) is not type(self.sourceImage):
                raise TypeError
            if np.ndim(containerImage) is not np.ndim(self.sourceImage):
                raise ValueError

            self.homography = Homography(sourcePoints=self.sourcePoints, targetPoints=self.targetPoints, effect=self.effect)
            x = [f[0] for f in self.targetPoints]
            y = [s[1] for s in self.targetPoints]

            rowmin, rowmax = min(y), max(y) + 1
            colmin, colmax = min(x), max(x) + 1



            if self.effect != None:
                self.sourcePoints = [(0,0), (self.sbound[1], 0), (0, self.sbound[0]), (self.sbound[1], self.sbound[0])]
                self.homography = Homography(sourcePoints=self.sourcePoints, targetPoints=self.targetPoints, effect=self.effect)
            else:
                self.homography = Homography(sourcePoints=self.sourcePoints, targetPoints=self.targetPoints, effect=self.effect)

            rbs =  interpolate.RectBivariateSpline(np.arange(0.0, self.sbound[0]), np.arange(0.0, self.sbound[1]), self.sourceImage, kx=1, ky=1)
            # Identify the target bounding box.
            # For every point within the box, perform an inverse projection of the coordinates.
            for i in np.arange(rowmin, rowmax):
                for j in np.arange(colmin, colmax):
                    a = j, i
                    px, py = self.homography.inverseProject(a)
                    if 0 <= px and px <= self.sbound[1] - 1 and 0 <= py and py <= self.sbound[0] - 1:
                        # read value of source and use 2D interpolation
                        containerImage[i][j] = rbs(py, px)

            return containerImage

def outer(A, B, ndim=1):
    """
    Computes outer product over the last axes of A and B.

    The other axes are broadcasted. Thus, if A has shape (..., N) and B has
    shape (..., M), then the result has shape (..., N, M).

    Using the argument `ndim` it is possible to change that how many axes
    trailing axes are used for the outer product. For instance, if ndim=3, A and
    B have shapes (...,N1,N2,N3) and (...,M1,M2,M3), the result has shape
    (...,N1,M1,N2,M2,N3,M3).
    """
    if not utils.is_integer(ndim) or ndim < 0:
        raise ValueError('ndim must be non-negative integer')
    if ndim > 0:
        if ndim > np.ndim(A):
            raise ValueError('Argument ndim larger than ndim of the first '
                             'parameter')
        if ndim > np.ndim(B):
            raise ValueError('Argument ndim larger than ndim of the second '
                             'parameter')
        shape_A = np.shape(A) + (1,)*ndim
        shape_B = np.shape(B)[:-ndim] + (1,)*ndim + np.shape(B)[-ndim:]
        A = np.reshape(A, shape_A)
        B = np.reshape(B, shape_B)
    return A * B

    def convert_to_firing_rate(self, calciumtraces):

        firingrate = []
        denoisedcalicumtrace = []
        for ii in xrange(0, np.size(calciumtraces, 1)):
            fluor = calciumtraces[:, ii]
            print 'Deconvolving..', ii

            deconvolvedresult = constrained_foopsi(fluor)
            print np.ndim(deconvolvedresult[5])

            if np.ndim(deconvolvedresult[5]) > 1:
                print 'Skipping..Cell ', ii
                continue
            firingrate.append(deconvolvedresult[5])
            denoisedcalicumtrace.append(deconvolvedresult[0])

        firingrate = np.asarray(np.vstack(firingrate)).T
        denoisedcalicumtrace = np.asarray(np.vstack(denoisedcalicumtrace)).T

        np.savetxt(os.path.join(self.WorkingDirectory, key + "_denoiseddata.csv"), denoisedcalicumtrace, delimiter=",")
        np.savetxt(os.path.join(self.WorkingDirectory, key + "_firingrate.csv"), firingrate, delimiter=",")

        plt.plot(np.clip(firingrate, 0, np.max(firingrate)), '.')
        plt.savefig('/Users/seetha/Desktop/test1.png')

 def Zo(s):
     if np.ndim(h) == 1:
         s = s[:, np.newaxis]
     if np.ndim(h) == 2:
         s = s[:, np.newaxis, np.newaxis]
                 
     return (hc * s + C(s) * h) / (hc + h)