Python sparse.bmat方法代码示例

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def _compute_jacobian(self, J_eq, J_ineq, s):
        if self.n_ineq == 0:
            return J_eq
        else:
            if spc.issparse(J_eq) or spc.issparse(J_ineq):
                # It is expected that J_eq and J_ineq
                # are already `csr_matrix` because of
                # the way ``BoxConstraint``, ``NonlinearConstraint``
                # and ``LinearConstraint`` are defined.
                J_eq = spc.csr_matrix(J_eq)
                J_ineq = spc.csr_matrix(J_ineq)
                return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
            else:
                S = np.diag(s)
                zeros = np.zeros((self.n_eq, self.n_ineq))
                # Convert to matrix
                if spc.issparse(J_ineq):
                    J_ineq = J_ineq.toarray()
                if spc.issparse(J_eq):
                    J_eq = J_eq.toarray()
                # Concatenate matrices
                return np.asarray(np.bmat([[J_eq, zeros],
                                           [J_ineq, S]])) 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def _compute_jacobian(self, J_eq, J_ineq, s):
        if self.n_ineq == 0:
            return J_eq
        else:
            if sps.issparse(J_eq) or sps.issparse(J_ineq):
                # It is expected that J_eq and J_ineq
                # are already `csr_matrix` because of
                # the way ``BoxConstraint``, ``NonlinearConstraint``
                # and ``LinearConstraint`` are defined.
                J_eq = sps.csr_matrix(J_eq)
                J_ineq = sps.csr_matrix(J_ineq)
                return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
            else:
                S = np.diag(s)
                zeros = np.zeros((self.n_eq, self.n_ineq))
                # Convert to matrix
                if sps.issparse(J_ineq):
                    J_ineq = J_ineq.toarray()
                if sps.issparse(J_eq):
                    J_eq = J_eq.toarray()
                # Concatenate matrices
                return np.asarray(np.bmat([[J_eq, zeros],
                                           [J_ineq, S]])) 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def _merge_blocks(blocks):
        """
        Helper function that merges the _blocks attribute of a ds-array into
        a single ndarray / sparse matrix.
        """
        sparse = None
        b0 = blocks[0][0]
        if sparse is None:
            sparse = issparse(b0)

        if sparse:
            ret = sp.bmat(blocks, format=b0.getformat(), dtype=b0.dtype)
        else:
            ret = np.block(blocks)

        return ret 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def edge_dos(intra0,inter0,scale=4.,w=20,npol=300,ne=500,bulk=False,
                use_random=True,nrand=20):
  """Calculated the edge DOS using the KPM"""
  h = [[None for j in range(w)] for i in range(w)]
  intra = csc_matrix(intra0)
  inter = csc_matrix(inter0)
  for i in range(w): h[i][i] = intra
  for i in range(w-1): 
    h[i+1][i] = inter.H
    h[i][i+1] = inter
  h = bmat(h) # sparse hamiltonian
  ds = np.zeros(ne)
  dsb = np.zeros(ne)
  norb = intra0.shape[0] # orbitals ina cell
  for i in range(norb):
    (xs,ys) = ldos0d(h,i=i,scale=scale,npol=npol,ne=ne) 
    ds += ys # store
    if bulk:
      (xs,zs) = ldos0d(h,i=w*norb//2 + i,scale=scale,npol=npol,ne=ne) 
      dsb += zs # store
  if not bulk: return (xs,ds/w)
  else: return (xs,ds/w,dsb/w) 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def add_interlayer(t,ri,rj,has_spin=True,is_sparse=True):
  """Calculate interlayer coupling"""
  m = np.matrix([[0. for i in ri] for j in rj])
  if has_spin: m = bmat([[csc(m),None],[None,csc(m)]]).todense()
  zi = [r[2] for r in ri]
  zj = [r[2] for r in rj]
  for i in range(len(ri)): # add the interlayer
    for j in range(len(ri)): # add the interlayer
      rij = ri[i] - rj[j] # distance between atoms
      dz = zi[i] - zj[j] # vertical distance
      if (3.99<rij.dot(rij)<4.01) and (3.99<(dz*dz)<4.01): # check if connect
        if has_spin: # spin polarized
          m[2*i,2*j] = t
          m[2*i+1,2*j+1] = t
        else:  # spin unpolarized
          m[i,j] = t
  return m 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def bulk2ribbon_zz(h,n=10):
  """Converts a hexagonal bulk hamiltonian into a ribbon hamiltonian"""
  h = honeycomb2square(h) # create a square 2d geometry
  go = h.geometry.copy() # copy geometry
  ho = h.copy() # copy hamiltonian
  ho.dimensionality = 1 # reduce dimensionality
  go.dimensionality = 1 # reduce dimensionality
  intra = [[None for i in range(n)] for j in range(n)]
  inter = [[None for i in range(n)] for j in range(n)]
  for i in range(n): # to the the sam index
    intra[i][i] = csc(h.intra) 
    inter[i][i] = csc(h.tx) 
  for i in range(n-1): # one more or less
    intra[i][i+1] = csc(h.ty)  
    intra[i+1][i] = csc(h.ty.H)  
    inter[i+1][i] = csc(h.txmy) 
    inter[i][i+1] = csc(h.txy) 
  ho.intra = bmat(intra).todense()
  ho.inter = bmat(inter).todense()
  return ho 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def ldos_finite(h,e=0.0,n=10,nwf=4,delta=0.0001):
  """Calculate the density of states for a finite system"""
  if h.dimensionality!=1: raise # if it is not one dimensional
  intra = csc(h.intra) # convert to sparse
  inter = csc(h.inter) # convert to sparse
  interH = inter.H # hermitian
  m = [[None for i in range(n)] for j in range(n)] # full matrix
  for i in range(n): # add intracell
    m[i][i] = intra
  for i in range(n-1): # add intercell
    m[i][i+1] = inter
    m[i+1][i] = interH
  m = bmat(m) # convert to matrix
  (ene,wfs) = slg.eigsh(m,k=nwf,which="LM",sigma=0.0) # diagonalize
  wfs = wfs.transpose() # transpose wavefunctions
  dos = (wfs[0].real)*0.0 # calculate dos
  for (ie,f) in zip(ene,wfs): # loop over waves
    c = 1./(1.+((ie-e)/delta)**2) # calculate coefficient
    dos += np.abs(f)*c # add contribution
  odos = spatial_dos(h,dos) # get the spatial distribution
  go = h.geometry.supercell(n) # get the supercell
  write_ldos(go.x,go.y,odos) # write in a file
  return dos # return the dos 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def block_inverse(m,i=0,j=0):
  """ Calculate a certain element of the inverse of a block matrix"""
  from scipy.sparse import csc_matrix,bmat
  nb = len(m) # number of blocks
  if i<0: i += nb 
  if j<0: j += nb 
  mt = [[None for ii in range(nb)] for jj in range(nb)]
  for ii in range(nb): # diagonal part
    mt[ii][ii] = csc_matrix(m[ii][ii])
  for ii in range(nb-1):
    mt[ii][ii+1] = csc_matrix(m[ii][ii+1])
    mt[ii+1][ii] = csc_matrix(m[ii+1][ii])
  mt = bmat(mt).todense() # create dense matrix
  # select which elements you need
  ilist = [m[ii][ii].shape[0] for ii in range(i)] 
  jlist = [m[jj][jj].shape[1] for jj in range(j)] 
  imin = sum(ilist)
  jmin = sum(jlist)
  mt = mt.I # calculate inverse
  imax = imin + m[i][i].shape[0]
  jmax = jmin + m[j][j].shape[1]
  mo = [ [mt[ii,jj] for jj in range(jmin,jmax)] for ii in range(imin,imax) ] 
  mo = np.matrix(mo)
  return mo 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def interface(h1,h2,k=[0.0,0.,0.],energy=0.0,delta=0.01):
  """Get the Green function of an interface"""
  from scipy.sparse import csc_matrix as csc
  from scipy.sparse import bmat
  gs1,sf1 = green_kchain(h1,k=k,energy=energy,delta=delta,
                   only_bulk=False,reverse=True) # surface green function 
  gs2,sf2 = green_kchain(h2,k=k,energy=energy,delta=delta,
                   only_bulk=False,reverse=False) # surface green function 
  #############
  ## 1  C  2 ##
  #############
  # Now apply the Dyson equation
  (ons1,hop1) = get1dhamiltonian(h1,k,reverse=True) # get 1D Hamiltonian
  (ons2,hop2) = get1dhamiltonian(h2,k,reverse=False) # get 1D Hamiltonian
  havg = (hop1.H + hop2)/2. # average hopping
  ons = bmat([[csc(ons1),csc(havg)],[csc(havg.H),csc(ons2)]]) # onsite
  self2 = bmat([[csc(ons1)*0.0,None],[None,csc(hop2@sf2@hop2.H)]])
  self1 = bmat([[csc(hop1@sf1@hop1.H),None],[None,csc(ons2)*0.0]])
  # Dyson equation
  ez = (energy+1j*delta)*np.identity(ons1.shape[0]+ons2.shape[0]) # energy
  ginter = (ez - ons - self1 - self2).I # Green function
  # now return everything, first, second and hybrid
  return (gs1,sf1,gs2,sf2,ginter) 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def bulk_and_surface(h1,nk=100,energies=np.linspace(-1.,1.,100),**kwargs):
  """Get the surface DOS of an interface"""
  from scipy.sparse import csc_matrix,bmat
  if h1.dimensionality==2:
      kpath = [[k,0.,0.] for k in np.linspace(0.,1.,nk)]
  else: raise
  ik = 0
  h1 = h1.get_multicell() # multicell Hamiltonian
  tr = timing.Testimator("DOS") # generate object
  dos_bulk = energies*0.0
  dos_sur = energies*0.0
  for k in kpath:
    tr.remaining(ik,len(kpath)) # generate object
    ik += 1
    outs = green.surface_multienergy(h1,k=k,energies=energies,**kwargs)
    dos_bulk += np.array([-algebra.trace(g[1]).imag for g in outs])
    dos_sur += np.array([-algebra.trace(g[0]).imag for g in outs])
  dos_bulk /= len(kpath)
  dos_sur /= len(kpath)
  np.savetxt("DOS.OUT",np.array([energies,dos_bulk,dos_sur]).T)
  return energies,dos_bulk,dos_sur 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def get_interface(h,fun=None):
  """Return an operator that projects onte the interface"""
  dind = 1 # index to which divide the positions
  if h.has_spin:  dind *= 2 # duplicate for spin
  if h.has_eh:  dind *= 2  # duplicate for eh
  iden = csc(np.matrix(np.identity(dind,dtype=np.complex))) # identity matrix
  r = h.geometry.r # positions
  out = [[None for ri in r] for rj in r] # initialize
  if fun is None: # no input function
    cut = 2.0 # cutoff
    if h.dimensionality==1: index = 1
    elif h.dimensionality==2: index = 2
    else: raise
    def fun(ri): # define the function
      if np.abs(ri[index])<cut: return 1.0
      else: return 0.0
  for i in range(len(r)): # loop over positions
    out[i][i] = fun(r[i])*iden 
  return bmat(out) # return matrix 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def get_si(h,i=1):
  """Return a certain Pauli matrix for the full Hamiltonian"""
  if not h.has_spin: return None # no spin
  if i==1: si = sx # sx matrix
  elif i==2: si = sy # sy matrix
  elif i==3: si = sz # sz matrix
  else: raise # unknown pauli matrix
  if h.has_eh: ndim = h.intra.shape[0]//4 # half the dimension
  else: ndim = h.intra.shape[0]//2 # dimension
  if h.has_spin: # spinful system
    op = [[None for i in range(ndim)] for j in range(ndim)] # initialize
    for i in range(ndim): op[i][i] = si # store matrix
    op = bmat(op) # create matrix
  if h.has_eh: op = build_eh(op,is_sparse=True) # add electron and hole parts 
  return op

# define the functions for the three spin components 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def eh_operator(m):
  """Return the electron hole symmetry operator, as a function"""
  n = m.shape[0]//4 # number of sites 
  from .hamiltonians import sy
  msy = [[None for ri in range(n)] for j in range(n)]
  for i in range(n): msy[i][i] = sy # add sy
  msy = bmat(msy) # sy matrix in the electron subspace
  out = [[None,1j*msy],[-1j*msy,None]]
#  out = [[None for i in range(n)] for j in range(n)] # output matrix
#  tau = csc_matrix([[0.,0.,0.,-1.],[0,0,1,0],[0,1,0,0],[-1,0,0,0]])
#  for i in range(n):
#    out[i][i] = tau
  out = bmat(out) # sparse matrix
  out = reorder(out) # reshufle the matrix
  def ehop(inm):
    """Function that applies electron-hole symmetry to a matrix"""
    return out*np.conjugate(inm)*out # return matrix
  return ehop # return operator 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def rotate_mag(m,angle=0.0,axis=[0,0,1]):
  """ Rotates the magnetization along a particular axis """
  nat = len(m)/2

  #####################
  ## rotate the dmat ##
  #####################
  from scipy.sparse import coo_matrix, bmat
  sx = np.matrix([[0.,1.],[1.,0.]])
  sy = np.matrix([[0.,1j],[-1j,0.]])
  sz = np.matrix([[1.,0.],[0.,-1.]])

  from scipy.linalg import expm
  rmat = expm(1j*np.pi*angle*(axis[0]*sx + axis[1]*sy + axis[2]*sz))
  rmat = coo_matrix(rmat) # to sparse
  rot = [[None for i in range(nat)] for j in range(nat)]
  for i in range(nat):
    rot[i][i] = rmat
  rot = bmat(rot).todense() # rotational matrix

  # rotate the matrix
  m = rot * m * rot.H

  return m 

# 需要导入模块: from scipy import sparse [as 别名]
# 或者: from scipy.sparse import bmat [as 别名]
def homogeneous_form(self):
        return sp.bmat([[self.P, self.q/2], [self.q.T/2, self.r]])

    # Returns QuadraticFunction f1, f2 such that
    # f(x) = f1(x) - f2(x), with f1 and f2 both convex.
    # Affine and constant components are always put into f1.