numpy基础属性方法随机整理(9):专用函数-- np.lexsort() / np.sort_complex两种方法实现间接联合排序

本文介绍了一种名为间接联合排序的方法,该方法通过获取排序样本的下标而非直接排序元素本身来实现联合排序。文章详细解释了如何利用NumPy库中的lexsort函数进行联合排序,并提供了一个使用复数进行联合排序的替代方案。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

间接联合排序: 间接获取排序样本的下标
原始数列:8 2 3 1 7 4 6 5 9
直接排序:1 2 3 4 5 6 7 8 9
间接排序:3 1 2 5 7 6 4 0 8 (原始序列元素的下标)
姓名:张三 李四 王五 赵六 陈七
成绩:90 70 50 80 60
下标:0 1 2 3 4
成绩升序对应的下标:2 4 1 3 0
年龄:20 30 30 20 40 (有相同值,对年龄排序时参考成绩进行联合排序)
下标:
联合排序对应的下标:3 0 2 1 4(年龄升序)

函数语法:

   numpy.lexsort((参考序列,待排序列)) ---> return:索引序列
   np.sort_complex(复数数组) ---> 按实部的升序排列,实部相同的参考虚部的升序排列
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 27 16:36:45 2018

@author: Administrator
"""

'''
np.where(namesComplexSorted == 'Kity')
    返回值:(array([1], dtype=int64),)   tuple和list复合的
    需要处理后才能传入list中用于后续使用 : val_return[0][0]
np.take(ages, np.where(ages == v))
    参数:(数组a, 下标数组)
    返回值:数组a中由下标数组所表示的a数组中的元素集合
'''

import numpy as np

names = np.array(['Tom','Mike','Jue','Kity','scote'])
scores = np.array([90., 70., 68., 80., 76.])
ages = np.array([29, 23, 20, 22, 27])               # [29 23 20 22 27]

# 联合排序方法1:返回下标index
index_lexsorted = np.lexsort((scores, ages))        # [2 3 1 4 0]
print(index_lexsorted)

names_lexsorted = names[np.lexsort((scores, ages))]
print(names_lexsorted)      # ['Jue' 'Kity' 'Mike' 'scote' 'Tom']

print('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')

# 联合排序方法2: 巧用复数进行联合排序
c = ages + scores * 1j                              # j要用 1j 表示
d = np.sort_complex(c).real                         # float:[20. 22. 23. 27. 29.] <class 'numpy.ndarray'>
d = d.astype('int64')                               # int:[20 22 23 27 29] <class 'numpy.ndarray'>

index_names = []
for i, v in enumerate(d):
    print(np.where(ages ==v)[0][0], np.take(ages, np.where(ages == v))[0][0])
    index_names.append(np.where(ages ==v)[0][0])    # 按照d的元素顺序依次对比,取出ages的下标

namesComplexSorted = names[index_names]
print(namesComplexSorted)                           # ['Jue' 'Kity' 'Mike' 'scote' 'Tom']
print(np.where(namesComplexSorted == 'Kity')[0][0]) # (array([1], dtype=int64),)

import numpy as np import matplotlib.pyplot as plt import matplotlib.cm as cm import matplotlib.colors as mcolors import matplotlib as mpl # 确保所有必需的函数都被定义 # --- 哈密顿量及其导数计算函数 --- def ham(lat,orb,hop,k): """ 计算给定k点的哈密顿量H(k)。 """ h=np.zeros((len(orb),len(orb)),dtype=&#39;complex&#39;) lat_np = np.array(lat) # 确保lat是numpy数组以便进行点积 for i in range (len(hop)): t_val = hop[i][0] # 跳跃项的值 from_orb_idx = hop[i][1] # 起始轨道索引 to_orb_idx = hop[i][2] # 结束轨道索引 Rvec_coords = np.array(hop[i][3:5]) # 晶格位移矢量 (rx, ry) # 计算真实空间中的位移矢量 R + r_to - r_from # (Rvec_x, Rvec_y) + (orb_to_x - orb_from_x, orb_to_y - orb_from_y) orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx]) latr = np.dot((Rvec_coords + orbital_diff), lat_np) # 转换为笛卡尔坐标系下的真实空间位移 newterm = t_val * np.exp(1j * (np.dot(latr, k))) h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm # 哈密顿量是厄米的,所以添加共轭转置部分 h = h + np.conjugate(h).transpose() return h def hamx(lat,orb,hop,k): """ 计算哈密顿量对kx的导数 dH/dkx。 """ h=np.zeros((len(orb),len(orb)),dtype=&#39;complex&#39;) lat_np = np.array(lat) for i in range (len(hop)): t_val = hop[i][0] from_orb_idx = hop[i][1] to_orb_idx = hop[i][2] Rvec_coords = np.array(hop[i][3:5]) orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx]) latr = np.dot((Rvec_coords + orbital_diff), lat_np) newterm = t_val * np.exp(1j * (np.dot(latr, k))) * (1j * latr[0]) # 乘以 i * (R_x + r_y_x - r_x_x) h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm h = h + np.conjugate(h).transpose() return h def hamy(lat,orb,hop,k): """ 计算哈密顿量对ky的导数 dH/dky。 """ h=np.zeros((len(orb),len(orb)),dtype=&#39;complex&#39;) lat_np = np.array(lat) for i in range (len(hop)): t_val = hop[i][0] from_orb_idx = hop[i][1] to_orb_idx = hop[i][2] Rvec_coords = np.array(hop[i][3:5]) orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx]) latr = np.dot((Rvec_coords + orbital_diff), lat_np) newterm = t_val * np.exp(1j * (np.dot(latr, k))) * (1j * latr[1]) # 乘以 i * (R_y + r_y_y - r_x_y) h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm h = h + np.conjugate(h).transpose() return h # --- 用户提供的模型参数 --- lat=np.array([[1.0,0.0],[0.0,1.0]]) orb=np.array([[0.0,0.0],[0.5,0.0],[0.0,0.5]]) t1=1.0 t2=0.4*np.exp((1.j)*np.pi*0.37) t3=0.13 # 将pythtb的set_hop定义转换为我们ham函数所需的hop列表格式 # 格式: [hoppping_term, orbital_from, orbital_to, lattice_vector_x, lattice_vector_y] hop=[ [t1, 0, 1, 0, 0], [t1, 0, 2, 0, 0], [t1, 0, 1, -1, 0], [t1, 0, 2, 0, -1], [t2, 1, 2, 0, 0], [t2, 2, 1, -1, 0], [t2, 1, 2, 1, -1], [t2, 2, 1, 0, 1], [t3, 0, 0, 1, 0], [t3, 0, 0, 0, 1], [t3, 1, 1, 1, 0], [t3, 1, 1, 0, 1], [t3, 2, 2, 1, 0], [t3, 2, 2, 0, 1], ] # 在位能 (onsite energies),根据原始代码设置为零 onsite=[0]*len(orb) # --- K空间网格设置 --- nx=200 # x方向K点数 ny=230 # y方向K点数 kxlim=2*np.pi # Kx方向的半范围 (-kxlim到+kxlim) kylim=2*np.pi # Ky方向的半范围 (-kylim到+kylim) # K空间总面积为 (2*kxlim) * (2*kylim) # ds 为每个K点代表的K空间面积元 ds=(4*kxlim*kylim)/(nx*ny) origin=[0,0] # K空间原点 # 初始化存储计算结果的数组 eccent=np.zeros((nx+1,ny+1)) # 特征椭球离心率 area=np.zeros((nx+1,ny+1)) # 特征椭球面积 curva=np.zeros((nx+1,ny+1)) # 贝里曲率 (未积分) differ=np.zeros((nx+1,ny+1)) # tr(g)-|B| kxlist=np.zeros(nx+1) kylist=np.zeros(ny+1) # --- K空间循环与计算 --- print("🚀 开始计算K空间数据...") for x in range(nx+1): kx=(x-nx/2)/(nx/2)*kxlim+origin[0] # 线性映射kx到 [-kxlim, kxlim] kxlist[x]=kx for y in range(ny+1): ky=(y-ny/2)/(ny/2)*kylim+origin[1] # 线性映射ky到 [-kylim, kylim] kylist[y]=ky # 1. 计算哈密顿量H(k) H=ham(lat,orb,hop,[kx,ky]) for i in range(len(H)): H[i,i]= H[i,i]+onsite[i] # 添加在位能 # 2. 计算哈密顿量对k的导数Hx和Hy Hx=hamx(lat,orb,hop,[kx,ky]) Hy=hamy(lat,orb,hop,[kx,ky]) # 3. 对角化H(k)得到本征值e和本征矢量v e,v=np.linalg.eigh(H) # e: eigenvalues, v: eigenvectors (columns are eigenvectors) # 4. 计算量子度规张量Q Q=np.zeros((2,2),dtype=&#39;complex&#39;) # 根据量子度规张量的定义 Q_ab = sum_{n!=m} <psi_m|d_a H|psi_n><psi_n|d_b H|psi_m> / (E_m - E_n)^2 # 这里计算的是基态 (m=0) 的量子度规张量 for i in range(1,len(orb)): # 遍历所有激发态 (n=1, 2, ...) # Q_xx Q[0,0]=Q[0,0]+np.dot(v[:,0].conjugate(),np.dot(Hx,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hx,v[:,0]))/((e[0]-e[i])**2) # Q_xy Q[0,1]=Q[0,1]+np.dot(v[:,0].conjugate(),np.dot(Hx,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hy,v[:,0]))/((e[0]-e[i])**2) # Q_yx Q[1,0]=Q[1,0]+np.dot(v[:,0].conjugate(),np.dot(Hy,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hx,v[:,0]))/((e[0]-e[i])**2) # Q_yy Q[1,1]=Q[1,1]+np.dot(v[:,0].conjugate(),np.dot(Hy,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hy,v[:,0]))/((e[0]-e[i])**2) # 5. 从Q中提取量子度规张量g和贝里曲率B g=Q.real # g 是量子度规张量的实部 (g_ab = Re(Q_ab)) F=-2*Q.imag # 贝里曲率张量 (F_ab = -2 Im(Q_ab)) B=F[0][1] # 贝里曲率分量 (F_xy) # 6. 计算特征椭球的参数 evals,evec=np.linalg.eigh(g) # g 的本征值 (代表半轴长的平方) c=np.sqrt(np.abs(evals[0]-evals[1])) # 焦点距离 a=np.sqrt(np.max(evals)) # 最长半轴 e=c/a if a != 0 else 0.0 # 离心率,避免除以零 # 7. 存储结果 eccent[x,y]=e area[x,y]=np.sqrt(np.linalg.det(g))*np.pi # 椭球面积 = pi * a * b = pi * sqrt(g1*g2) = pi * sqrt(det(g)) curva[x,y]=B*ds # 贝里曲率密度 differ[x,y]=(np.trace(g)-np.abs(B))*ds # tr(g)-|B| 密度 print("✅ K空间数据计算完成!") # --- 绘制图像 --- print("📊 正在生成图表...") fig, ([ax1,ax2],[ax3,ax4]) = plt.subplots(2,2,figsize=(12,12)) # 1. 贝里曲率图 (Berry Curvature) ax1.imshow(curva.T, origin="lower", cmap=cm.RdBu_r, extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()]) # 计算 Chern Number (假设K空间范围是4个布里渊区) chen = np.sum(curva) / (2 * np.pi) / 4 norm = mpl.colors.Normalize(vmin=curva.min(), vmax=curva.max()) cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r) # 使用RdBu_r以使正值显示红色,负值显示蓝色 ax1.set_title(f&#39;Berry Curvature (C = {round(chen, 3)})&#39;) # 标题显示切恩数 fig.colorbar(cmap, ax=ax1, orientation=&#39;vertical&#39;, shrink=0.945) ax1.set_xlabel(&#39;$k_x$&#39;) ax1.set_ylabel(&#39;$k_y$&#39;) ax1.set_aspect(&#39;equal&#39;, adjustable=&#39;box&#39;) # 保持宽高比 # 2. tr(g)-|B| 图 ax2.imshow(differ.T, origin="lower", cmap=cm.RdBu_r, extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()]) norm = mpl.colors.Normalize(vmin=differ.min(), vmax=differ.max()) cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r) ax2.set_title(&#39;$\\mathrm{tr}(g)-|B|$&#39;) # LaTeX格式的标题 fig.colorbar(cmap, ax=ax2, orientation=&#39;vertical&#39;, shrink=0.945) ax2.set_xlabel(&#39;$k_x$&#39;) ax2.set_ylabel(&#39;$k_y$&#39;) ax2.set_aspect(&#39;equal&#39;, adjustable=&#39;box&#39;) # 3. 特征椭球面积图 (Characteristic ellipsoid area) ax3.imshow(area.T, origin="lower", cmap=cm.RdBu_r, extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()]) norm = mpl.colors.Normalize(vmin=area.min(), vmax=area.max()) cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r) ax3.set_title(&#39;Characteristic Ellipsoid Area&#39;) fig.colorbar(cmap, ax=ax3, orientation=&#39;vertical&#39;, shrink=0.945) ax3.set_xlabel(&#39;$k_x$&#39;) ax3.set_ylabel(&#39;$k_y$&#39;) ax3.set_aspect(&#39;equal&#39;, adjustable=&#39;box&#39;) # 4. 特征椭球离心率图 (Characteristic ellipsoid eccentricity index) ax4.imshow(eccent.T, origin="lower", cmap=cm.RdBu_r, extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()]) norm = mpl.colors.Normalize(vmin=eccent.min(), vmax=eccent.max()) cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r) ax4.set_title(&#39;Characteristic Ellipsoid Eccentricity Index&#39;) fig.colorbar(cmap, ax=ax4, orientation=&#39;vertical&#39;, shrink=0.945) ax4.set_xlabel(&#39;$k_x$&#39;) ax4.set_ylabel(&#39;$k_y$&#39;) ax4.set_aspect(&#39;equal&#39;, adjustable=&#39;box&#39;) plt.tight_layout() # 自动调整子图参数,使之填充整个图像区域 plt.show() print("🎉 图表已成功显示!") 这有什么问题吗
最新发布
07-11
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值