def _expandg(g):
    """ Expand transition function to a matrix
    """
    P = sparse.coo_matrix((sp.ones(sp.prod(g.shape)),
                             (sp.r_[0:sp.prod(g.shape)],
                              g.flatten(1))))
    return P.tocsr()

	def step(self, x, y):
		SIZE = self.SIZE
		J = self.J
		h = self.h
		
			
		
		factor = y%2 * 2 - 1
		neighbours = scipy.array([
						self.m[x][(y + 1)% SIZE ] , 
						self.m[x][(y - 1)% SIZE ] ,
						self.m[(x+1)%SIZE][y] ,
						self.m[(x-1)%SIZE][y] ,
						self.m[(x-factor)%SIZE][(y-1)%SIZE] ,
						self.m[(x-factor)%SIZE][(y+1)%SIZE]
					    ])
			    
		jump = random.randint(1,int(35*self.cutoff)+1)*5 * random.choice([1,-1])
		
		diff0 = (neighbours - self.m[x][y])%360
		diff1 = (neighbours - self.m[x][y] - jump)%360
		
		p =  scipy.prod([self.be[i] for i in diff1]) / scipy.prod([self.be[i] for i in diff0])
		
		self.sts += 1
		if random.random() < p:
			self.change += 1
			self.m[x][y] = (self.m[x][y] + jump)%360
			self.E = self.E + (sum([self.cos[i] for i in diff1]) * J * (-1) - sum([self.cos[i] for i in diff0]) * J * (-1))

def genLine(grid, numElements):
  for seg in genNegSlopeDiagonal(grid, numElements):
    yield (scipy.prod(seg), seg, 'negative slope')
  for seg in genPosSlopeDiagonal(grid, numElements):
    yield (scipy.prod(seg), seg, 'positive slope')
  for seg in genHoriz(grid, numElements):
    yield (scipy.prod(seg), seg, 'horizontal')
  for seg in genVert(grid, numElements):
    yield (scipy.prod(seg), seg, 'vertical')

def cube_grid(dims):
    """
    Return a regular nD-cube mesh with given shape.

    Eg.
      cube_grid_nd((2,2))   -> 2x2   - 2d mesh (x,y)
      cube_grid_nd((4,3,2)) -> 4x3x2 - 3d mesh (x,y,z)

    Eg.
    
      v,i = cube_grid_nd((2,1))

      v =
      array([[ 0.,  0.],
             [ 1.,  0.],
             [ 2.,  0.],
             [ 0.,  1.],
             [ 1.,  1.],
             [ 2.,  1.]])

      i = 
      array([[[0, 3],
              [1, 4]],

             [[1, 4],
              [2, 5]]])

    """
    dims = tuple(dims)
    
    vert_dims = tuple(x+1 for x in dims)
    N = len(dims)
    
    vertices = zeros((prod(vert_dims),N))
    grid     = mgrid[tuple(slice(0,x,None) for x in reversed(vert_dims))]
    for i in range(N):
        vertices[:,i] = ravel(grid[N-i-1])


    #construct one cube to be tiled
    cube  = zeros((2,)*N,dtype='i')
    cycle = array([1] + list(cumprod(vert_dims)[:-1]),dtype='i')
    for i in ndindex(*((2,)*N)):
        cube[i] = sum(array(i) * cycle)
        cycle = array([1] + list(cumprod(vert_dims)[:-1]),dtype='i')


    #indices of all vertices which are the lower corner of a cube
    interior_indices = arange(prod(vert_dims)).reshape(tuple(reversed(vert_dims))).T
    interior_indices = interior_indices[tuple(slice(0,x,None) for x in dims)]

    indices = tile(cube,(prod(dims),) + (1,)*N) + interior_indices.reshape((prod(dims),) + (1,)*N)
    
    return (vertices,indices)

    def createRules(self, cages):
        """
        Create a random set of rules for the latin square given by
        self.solution under the cage structure defined in cages
        """

        for cage in cages:

            if len(cage) == 1:
                operation = EQUALS
                value = self.solution[cage[0]]

            else:

                if len(cage) == 2:
                    operation = random.randint(1, 4)
                    values = sorted([self.solution[index] for index in cage])

                    if operation == PLUS:
                        value = scipy.sum(values)

                    elif operation == MINUS:
                        value = values[1] - values[0]

                    elif operation == TIMES:
                        value = scipy.prod(values)

                    elif operation == DIVIDE:
                        vMin, vMax = values
                        if vMax % vMin == 0:
                            value = vMax / vMin
                        else:
                            operation = MINUS
                            value = vMax - vMin

                    else:
                        raise ValueError, "Operation -- I'm the doctor for you"

                else:
                    operation = random.randint(1, 2)
                    values = sorted([self.solution[index] for index in cage])

                    if operation == PLUS:
                        value = scipy.sum(values)

                    elif operation == TIMES:
                        value = scipy.prod(values)

                    else:
                        raise ValueError, "Operation not possible for this list, dog"

            self.ruleList.append((cage, operation, value))

def diffmat(dims,order = 'C'):
    """ This function will return a tuple of difference matricies for data from an 
        Nd array that has been rasterized. The order parameter determines whether 
        the array was rasterized in a C style (python) of FORTRAN style (MATLAB).
        Inputs:
            dims- A list of the size of the x,y,z.. dimensions.
            order- Specifies the vectorization of the matrix
        Outputs:
            dx,dy,dy... - The finite difference operators for a vectorized array.
                If these are to be stacked together as one big operator then
                sp.sparse.vstack should be used.
    """
    # flip the dimensions around
    dims=[int(i) for i in dims]
    xdim = dims[0]
    ydim = dims[1]
    dims[0]=ydim
    dims[1]=xdim
    
    
    if order.lower() == 'c':
        dims = dims[::-1]

    outD = []
    for idimn, idim in enumerate(dims):
        if idim==0:
            outD.append(sp.array([]))
            continue
        e = sp.ones(idim)
        dthing = sp.vstack((-e,e))
        D = sp.sparse.spdiags(dthing,[0,1],idim-1,idim).toarray()
        D = sp.vstack((D,D[-1]))
        if idim>0:
            E = sp.sparse.eye(sp.prod(dims[:idimn]))
            D = sp.sparse.kron(D,E)

        if idimn<len(dims)-1:
            E = sp.sparse.eye(sp.prod(dims[idimn+1:]))
            D = sp.sparse.kron(E,D)

        outD.append(sp.sparse.csc_matrix(D))
    if order.lower() == 'c':
        outD=outD[::-1]
    Dy=outD[0]
    Dx = outD[1]
    outD[0]=Dx
    outD[1]=Dy
    return tuple(outD)

def _steadystate_eigen(L, ss_args):
    """
    Internal function for solving the steady state problem by
    finding the eigenvector corresponding to the zero eigenvalue
    of the Liouvillian using ARPACK.
    """
    if settings.debug:
        print('Starting Eigen solver...')

    dims = L.dims[0]
    shape = prod(dims[0])
    L = L.data.tocsc()

    if ss_args['use_rcm']:
        if settings.debug:
            old_band = sp_bandwidth(L)[0]
            print('Original bandwidth:', old_band)
        perm = reverse_cuthill_mckee(L)
        rev_perm = np.argsort(perm)
        L = sp_permute(L, perm, perm, 'csc')
        if settings.debug:
            rcm_band = sp_bandwidth(L)[0]
            print('RCM bandwidth:', rcm_band)
            print('Bandwidth reduction factor:', round(old_band/rcm_band, 1))

    eigval, eigvec = eigs(L, k=1, sigma=1e-15, tol=ss_args['tol'],
                          which='LM', maxiter=ss_args['maxiter'])

    if ss_args['use_rcm']:
        eigvec = eigvec[np.ix_(rev_perm,)]

    data = vec2mat(eigvec)
    data = 0.5 * (data + data.conj().T)
    out = Qobj(data, dims=dims, isherm=True)
    return out/out.tr()

    def _generate_pores(self):
        r"""
        Generate the pores (coordinates, numbering and types)
        """
        self._logger.info("generate_pores: Create specified number of pores")

        #Find non-zero elements in image
        template = self._template
        Np = np.sum(template > 0)
        #Add pores to data and ifo
        pind = np.arange(0, Np)
        self.set_pore_info(label='all', locations=pind)
        self.set_pore_data(prop='numbering', data=pind)  # Remove eventually

        
        img_ind = np.ravel_multi_index(sp.nonzero(template), dims=sp.shape(template), order='F')
        self.set_pore_data(prop='voxel_index', data=img_ind)

        #This voxel_to_pore map is messy but works
        temp = sp.prod(sp.shape(template))*sp.ones(np.prod(sp.shape(template),),dtype=sp.int32)
        temp[img_ind] = pind
        self._voxel_to_pore_map = temp

        coords = self._Lc*(0.5 + np.transpose(np.nonzero(template)))
        self.set_pore_data(prop='coords', data=coords)
        self._logger.debug("generate_pores: End of method")

    def fit(self):

        self.initialiseMetric()

        xMin, xMax, xStp, xShp = self.calcExhaustiveSearchGrid()
        numSearchPts = sp.prod(xShp)
    
        rootLogger.info("Grid search:")
        rootLogger.info("parameter x min   = %s" % (xMin,))
        rootLogger.info("parameter x max   = %s" % (xMax,))
        rootLogger.info("parameter x step  = %s" % (xStp,))
        rootLogger.info("parameter x shape = %s, %s metric evaluations" % (xShp, numSearchPts))
        rootLogger.info("Exhausive search...")
        cylList = []
        for cylIdx in range(0, self.numcyl):
            self.maskGradientImageCylinders(cylList)
            resultList = self.calcExhaustiveSearchResults(xMin, xMax, xStp, xShp)
            resultList = self.eliminatePoorResults(resultList)
            rootLogger.info("Done exhausive search.")

            resultList = self.calcBestRefinements(resultList)
            cylList.append(resultList[0])

        # Convert the parameter-vectors into 3 element centre-point, 3-element axis, etc.
        cylList = \
            [
                [resultPair[0], self.calcFullCylinderParameters(resultPair[1])]
                for
                resultPair in cylList
            ]

        return cylList

    def calcExhaustiveSearchResults(self, xMin, xMax, xStp, xShp):
        resultList = []
        numSearchPts = sp.prod(xShp)
        if (self.distributedMetricEvaluation):
            if ((self.comm == None) or (self.comm.Get_rank() == self.root)):
                for i in range(0, numSearchPts):
                    x = xMin + np.unravel_index(i, xShp)*xStp
                    resultList.append([self.metric(x), x])
                    rootLogger.debug(resultList[-1])
                self.metric.rootTerminate()
            else:
                self.metric.waitForEvaluate()
            if (self.comm != None):
                resultList = self.comm.bcast(resultList, self.root) 
        else:
            commSz = 1
            commRk = 0
            if (self.comm != None):
                commSz = self.comm.Get_size()
                commRk = self.comm.Get_rank()
            for i in range(commRk, numSearchPts, commSz):
                x = xMin + np.unravel_index(i, xShp)*xStp
                resultList.append([self.metric(x), x])
            if (self.comm != None):
                rListList = self.comm.allgather(resultList)
                resultList = []
                for rList in rListList:
                    resultList += rList

        return resultList

def state_number_index(dims, state):
    """
    Return the index of a quantum state corresponding to state,
    given a system with dimensions given by dims.

    Example:

        >>> state_number_index([2, 2, 2], [1, 1, 0])
        6

    Parameters
    ----------
    dims : list or array
        The quantum state dimensions array, as it would appear in a Qobj.

    state : list
        State number array.

    Returns
    -------
    idx : int
        The index of the state given by `state` in standard enumeration
        ordering.

    """
    return int(
        sum([state[i] * prod(dims[i + 1:]) for i, d in enumerate(dims)]))

def _steadystate_direct_sparse(L, verbose=False):
    """
    Direct solver that use scipy sparse matrices
    """
    if verbose:
        print('Starting direct solver...')

    n = prod(L.dims[0][0])
    b = sp.csr_matrix(([1.0], ([0], [0])), shape=(n ** 2, 1), dtype=complex)
    M = L.data + sp.csr_matrix((np.ones(n),
            (np.zeros(n), [nn * (n + 1) for nn in range(n)])),
            shape=(n ** 2, n ** 2))
    
    use_solver(assumeSortedIndices=True, useUmfpack=False)
    M.sort_indices()

    if verbose:
        start_time = time.time()
    # Do the actual solving here
    v = spsolve(M, b)

    if verbose:
        print('Direct solver time: ', time.time() - start_time)
    
    data = vec2mat(v)
    data = 0.5 * (data + data.conj().T)

    return Qobj(data, dims=L.dims[0], isherm=True)

def _steadystate_direct_dense(L, verbose=False):
    """
    Direct solver that use numpy dense matrices. Suitable for
    small system, with a few states.
    """
    if verbose:
        print('Starting direct dense solver...')

    n = prod(L.dims[0][0])
    b = np.zeros(n ** 2)
    b[0] = 1.0

    M = L.data.todense()
    M[0, :] = np.diag(np.ones(n)).reshape((1, n ** 2))
    if verbose:
        start_time = time.time()
    v = np.linalg.solve(M, b)

    if verbose:
        print('Direct dense solver time: ', time.time() - start_time)

    data = vec2mat(v)
    data = 0.5 * (data + data.conj().T)

    return Qobj(data, dims=L.dims[0], isherm=True)

def _steadystate_lu(L, use_rcm=True, use_umfpack=False):
    """
    Find the steady state(s) of an open quantum system by computing the
    LU decomposition of the underlying matrix.
    """
    if settings.debug:
        print('Starting LU solver...')
    dims=L.dims[0]
    weight=np.abs(L.data.max())
    n = prod(L.dims[0][0])
    b = np.zeros(n ** 2, dtype=complex)
    b[0] = weight
    L = L.data.tocsc() + sp.csc_matrix((weight*np.ones(n),
                    (np.zeros(n), [nn * (n + 1) for nn in range(n)])),
        shape=(n ** 2, n ** 2))
    
    L.sort_indices()
    use_solver(assumeSortedIndices=True, useUmfpack=use_umfpack)
    if use_rcm:
        perm = symrcm(L)
        L = sparse_permute(L,perm,perm)
        b = b[np.ix_(perm,)]
    
    solve = factorized(L)
    v = solve(b)
    if use_rcm:
        rev_perm = np.argsort(perm)
        v = v[np.ix_(rev_perm,)]
    data = vec2mat(v)
    data = 0.5 * (data + data.conj().T)

    return Qobj(data, dims=dims, isherm=True)

 def __voxel_4conectedness(self, shape):
     """
     Returns the number of edges for the supplied image shape assuming 4-connectedness.
     """
     shape = list(shape)
     while 1 in shape: shape.remove(1)
     return int(round(sum([(dim - 1)/float(dim) for dim in shape]) * scipy.prod(shape)))

def isket(Q):
    """
    Determines if given quantum object is a ket-vector.
	
	Parameters
	----------
	Q : qobj
	    Quantum object
	
	Returns
	------- 
	isket : bool
	    True if qobj is ket-vector, False otherwise.
	
	Examples
	--------	    
    >>> psi=basis(5,2)
    >>> isket(psi)
    True
	    
	"""
    result = isinstance(Q.dims[0],list)
    if result:
        result = result and prod(Q.dims[1])==1
    return result

def hist_flatened(im, nbr_bins=10):
    """
    @param im: the (gray-scale) image as numpy/scipy array
    @param nbr_bins: the number of bins
    @return: the bins of the flattened histogram of the image
    """
    # get image histogram
    imhist, bins = histogram(im.flatten(), 1000)

    # only take bins with content into account
    nz = imhist.nonzero()
    imhist = imhist[nz]
    bins = bins[nz]

    # prepare iteration
    bins_final = [bins[0]]  # set initial bin delimiter
    bins_content = scipy.prod(im.shape) / float(nbr_bins)
    tmp_content = 0
    for i in range(len(imhist) - 1):
        tmp_content += imhist[i]
        if tmp_content >= bins_content:  # bin full
            # bins_final.append(bins[i+1]) # add new bin delimiter
            # tmp_content = 0
            div = float(imhist[i]) / (bins_content - (tmp_content - imhist[i]))  # what i got / what i want
            bins_final.append(
                bins[i] + (bins[i + 1] - bins[i]) / div
            )  # append a partial bin border, assuming that the dist inside the bin in equal
            tmp_content = imhist[i] - (bins_content - (tmp_content - imhist[i]))

    bins_final.append(im.max() + 1)  # one added to work against rounding errors

    return bins_final

def _steadystate_direct_sparse(L, use_rcm=True, use_umfpack=False):
    """
    Direct solver that uses scipy sparse matrices
    """
    if settings.debug:
        print('Starting direct solver...')
    dims=L.dims[0]
    weight=np.abs(L.data.max())
    n = prod(L.dims[0][0])
    b = np.zeros((n ** 2, 1), dtype=complex)
    b[0,0] = weight
    L = L.data + sp.csr_matrix((weight*np.ones(n), (np.zeros(n), [nn * (n + 1) for nn in range(n)])),
                               shape=(n ** 2, n ** 2))
    L.sort_indices()
    use_solver(assumeSortedIndices=True, useUmfpack=use_umfpack)
    if use_rcm:
        perm = symrcm(L)
        L = sparse_permute(L,perm,perm)
        b = b[np.ix_(perm,)]
    
    v = spsolve(L, b)
    if use_rcm:
        rev_perm = np.argsort(perm)
        v = v[np.ix_(rev_perm,)]
    
    data = vec2mat(v)
    data = 0.5 * (data + data.conj().T)
    return Qobj(data, dims=dims, isherm=True)

 def compactRepresentation(self, photo):
     arr = scipy.ones((self.compareSize[0], self.compareSize[1], 3))*scipy.NaN
     if len(photo.shape) == 3: # an MxNx3 array of of ANY size:
         arr = photo
     elif photo.size == scipy.prod(self.fullSize)*3: # must be NxNx3, with N as in self.fullSize
         arr = photo.reshape((self.fullSize[0],self.fullSize[1],3))
     return scipy.misc.imresize(arr, self.compareSize).reshape((1,self.totalSize*3))

def steadystate_iterative(H, c_ops, use_precond=True):
    """
    .. note:: Experimental.
    """
    L = liouvillian_fast(H, c_ops)
    n = prod(L.dims[0][0])
    b = np.zeros(n ** 2)
    b[0] = 1.0
    A = L.data + sp.csr_matrix((np.ones(n), (np.zeros(n), \
            [nn * (n + 1) for nn in range(n)])), shape=(n ** 2, n ** 2))

    if use_precond:
        try:
            P = spilu(A, permc_spec='MMD_AT_PLUS_A')
            P_x = lambda x: P.solve(x)
            M = LinearOperator((n ** 2, n ** 2), matvec=P_x)
        except:
            warnings.warn("Preconditioning failed. Continuing without.",
                          UserWarning)
            M = None
    else:
        M = None

    v, check = bicgstab(A, b, tol=1e-5, M=M)

    return Qobj(vec2mat(v), dims=L.dims[0], isherm=True)