基于Numpy实现的仿真量子计算机（使用说明与源码）

Nervous_46216553

已于 2023-01-27 09:23:40 修改

阅读量668

点赞数

分类专栏：量子计算文章标签： numpy python

于 2023-01-26 23:12:07 首次发布

本文链接：https://blog.csdn.net/weixin_46216553/article/details/128767476

版权

量子计算专栏收录该内容

1 篇文章 0 订阅

订阅专栏

基于Numpy实现的仿真量子计算机

仿真器基本介绍
仿真器完整代码
演示代码
其它说明

仿真器基本介绍

可以保证10位机仿真实验的秒级响应速度。

数据类型介绍

量子态：np.array([[α],[β]])，α、β为复数但α的虚部是零，两者模长的平方之和与1的差在1e-16数量级
态空间：np.array([[a0],[a1],[a2],…,[an]]),a1到an都是复数，n=2^x-1,所有元素模长的平方和与1的差在1e-16数量级

量子位的定义函数介绍

zero_qubits(N) 获取由N个|0>态量子态组成的向量
one_qubits(N) 获取由N个|1>态量子态组成的向量
rand_qubits(N) 调试用，获取由N个随机态量子态组成的向量(α在(0,1)上均匀分布，β与实轴的角度在(0,2pi)上均匀分布，α=0时β=1)
inc_qubits(N) 调试用，(α从0递增到1，β与实轴的角度在(0,2pi)上均匀分布，α=0时β=1)

态空间获取函数

get_space(v) v是由多个量子态组成的向量，函数返回他们的张量积

操作门与操作函数介绍

1、提供基本单比特门：I，X，Y，Z，H，S，T
2、提供单比特旋转门：rx(t),ry(t),rz(t)和与rz等效但不会引入全局相位的相位门：ph(t)，t为角度
3、(不常用，内部函数用的)提供单比特控制门合成函数：C00(g),C01(g),C10(g),C11(g),CMN中M为控制位，N为控制方式,g为受控操作
C01(X)就是CNOT门，C01(C01(X))就是C^2NOT门
4、(不常用，内部函数用的)提供多比特门合成函数：G(g1,g2)，比如G(I,X)可以给两个量子位中的第二个作用X门
5、通用单bit门GN(i,n,g)：i作用在第几个bit，n总位数，g操作门
6、通用单bit控制门C1(i,j,n,g)：|1>有效 i控制bit，j受控bit，n总位数，g操作门，C0()作用相同|0>有效
注：i与j不可以相等，也不能越界
7、量子位交换操作：swap(i,j,space),space是态空间，函数直接返回交换后的态空间提高效率。
注：i与j不可以相等，也不能越界

实用工具介绍

unify(space)去除量子态与态空间的全局相位，便于数据分析，减少运算的bug
get_v(space):把大的态空间尽力拆成小的态空间，便于数据分析
注：返回由np数组组成的列表，用如下代码打印比较美观

for v in get_v(space):
    print(v)

get_unmusk_order(space):获取无遮挡的量子位顺序，是个整数列表，里面是从0开始的“下标”
如果第一位与第三位纠缠，则第二位即使是纯态也不会被get_v分离出来，所以有时需要去遮挡
reorder(space,order):根据order重新排列量子位并返回新的态空间
get_p(space)：各态出现的概率，最左边是状态[1,0,0,0,…,0]出现的概率，最右边是[0,0,0,0,…,1]出现的概率，比如量子态[[1],[0]]与[[1],[0]]张成状态[1,0,0,0]。

仿真器完整代码

import numpy as np
np.set_printoptions(precision=3, suppress=True)

#创建N个qubit
#置零的
def zero_qubits(N):
    return np.array([[[1+0j],[0j]]]*N)
#置一的（调试用）
def one_qubits(N):
    return np.array([[[0j],[1+0j]]]*N)
#随机的（调试用）
def rand_qubits(N):
    a=np.random.rand(N,1,1)
    b=np.sqrt(np.ones_like(a)-a**2)*np.exp(np.random.rand()*np.pi*2j)
    b[a==0]=1
    return np.concatenate((a,b),axis=1)
#“顺序增”的（调试用）
def inc_qubits(N):
    a=np.array(range(N))/N
    a=a.reshape(N,1,1)
    b=np.sqrt(np.ones_like(a)-a**2)*np.exp(np.random.rand()*np.pi*2j)
    b[a==0]=1
    return np.concatenate((a,b),axis=1)

#获取由多个态空间张成的态空间
def get_space(v):
    if(len(v)==1):
        return v[0]
    space=v[0]
    for q in v[1:]:
        space=np.kron(space,q)
    return space

#全局相位置零，便于观察比较
def unify(space):
    for i in range(len(space)):
        if abs(space[i,0])>0:
            space=space*(abs(space[i,0])/space[i,0])
            break
    return space

#获取由态空间的最小子空间组成的向量
def get_v(space):
    space=unify(space)
    curl=2 #当前待定子空间最小长度
    v=[]  #子空间向量
    while curl<len(space):
        succ=True
        curv=np.array([])
        gap=len(space)//curl
        for i in range(gap):
            absv=np.sqrt(np.dot(space[i::gap,0].conjugate().T,space[i::gap,0]).real)
            if absv>0:
                if len(curv)==0:
                    curv=unify(space[i::gap]/absv)
                if len(curv)>0 and np.dot(curv.conjugate().T,unify(space[i::gap]/absv))[0,0].real<1-1e-12:
                    succ=False
        if succ:
            maxi=list(abs(curv)).index(max(abs(curv)))
            space=unify(space[gap*maxi:gap*maxi+gap]/curv[maxi])
            v=v+[curv]
            curl=2
        else:
            curl=curl*2
            succ=True
    v=v+[space]
    return v

#获取态空间各测量值的概率
def get_p(space):
    return (abs(space)**2).T[0]

#基本单bit门
I = np.array([[1,0],
              [0,1]])
X = np.array([[0,1],
              [1,0]])
Y = np.array([[0,-1j],
              [1j,0 ]])
Z = np.array([[1,0 ],
              [0,-1]])
H = np.sqrt(2)/2*(Z+X)

S = np.array([[1,0 ],
              [0,1j]])
T = np.array([[1,0                 ],
              [0,np.exp(np.pi*1j/4)]])
#旋转门与相位门
rx = lambda theta : np.cos(theta/2)*I-1j*np.sin(theta/2)*X
ry = lambda theta : np.cos(theta/2)*I-1j*np.sin(theta/2)*Y
rz = lambda theta : np.cos(theta/2)*I-1j*np.sin(theta/2)*Z
ph = lambda theta : \
        np.array([[1,0               ],
                  [0,np.exp(theta*1j)]])

#两个计算基的外积
M00 = np.dot([[1],[0]],[[1,0]])
M11 = np.dot([[0],[1]],[[0,1]])

#CU，控制位，控制方式
C00  = lambda g:np.kron(M00,g)+np.kron(M11,np.eye(g.shape[0]))
C01  = lambda g:np.kron(M00,np.eye(g.shape[0]))+np.kron(M11,g)
C10  = lambda g:np.kron(g,M00)+np.kron(np.eye(g.shape[0]),M11)
C11  = lambda g:np.kron(np.eye(g.shape[0]),M00)+np.kron(g,M11)
#门的位数扩展与多比特门合成
G   = lambda g1,g2:np.kron(g1,g2)
#通用扩展门：作用位，总位数，操作
def GN(i,n,g):
    A=g
    for ii in range(i):
        A=G(I,A)
    for ii in range(i+1,n):
        A=G(A,I)
    return A
#通用|0>控制门：控制位，受控位，总位数，操作
def C0(i,j,n,g):
    if j>i:
        d=j-i
        A=GN(d-1,d,g)
        A=C00(A)
        A=GN(i,i+1,A)
        A=GN(0,n-j,A)
    if j<i: 
        i,j=j,i
        d=j-i
        A=GN(0,d,g)
        A=C10(A)
        A=GN(i,i+1,A)
        A=GN(0,n-j,A)
    return A
#通用|1>控制门：控制位，受控位，总位数，操作
def C1(i,j,n,g):
    if j>i:
        d=j-i
        A=GN(d-1,d,g)
        A=C01(A)
        A=GN(i,i+1,A)
        A=GN(0,n-j,A)
    if j<i: 
        i,j=j,i
        d=j-i
        A=GN(0,d,g)
        A=C11(A)
        A=GN(i,i+1,A)
        A=GN(0,n-j,A)
    return A
#通用SWAP操作：直接操作节省开销
def swap(i,j,space):
    N=int(np.log2(len(space)))
    A=C1(i,j,N,X)
    B=C1(j,i,N,X)
    space=A@space
    space=B@space
    space=A@space
    return space

#获取态空间的无遮挡的序列
def get_unmusk_order(space):
    bits_n = int(np.log2(len(space))) #位数
    graph  = np.zeros([bits_n]*2)     #图
    visit  = np.zeros([bits_n])       #访问标记
    order  = [0]*bits_n               #无遮挡顺序
    #################################################################
    #先看z测量独立性
    ps     = get_p(space)             #统计图
    pb     = [0]*bits_n               #每位是1态的概率
    index  = np.array(range(len(space)))
    for i in range(bits_n):           #获取每位是1态的概率
        pass
        select_value=2**i
        select_bits =(index&select_value)//select_value
        pb[i]       =sum(ps*select_bits)
    for i in range(bits_n):
        select_value_i=2**i
        select_bits_i =(index&select_value_i)//select_value_i
        graph[i,i]=1
        for j in range(i+1,bits_n):
            select_value_j=2**j
            select_bits_j =(index&select_value_j)//select_value_j
            select_bits   =select_bits_i&select_bits_j
            if abs(sum(ps*select_bits)-pb[i]*pb[j])>1e-12:
                graph[bits_n-i-1,bits_n-j-1]=1
                graph[bits_n-j-1,bits_n-i-1]=1
    #################################################################
    #再看x测量独立性
    g=ry(-np.pi/2)
    for i in range(1,bits_n):
        g=G(ry(-np.pi/2),g)
    ps     = get_p(g@space)           #统计图
    pb     = [0]*bits_n               #每位是1态的概率
    index  = np.array(range(len(space)))
    for i in range(bits_n):           #获取每位是1态的概率
        pass
        select_value=2**i
        select_bits =(index&select_value)//select_value
        pb[i]       =sum(ps*select_bits)
    for i in range(bits_n):
        select_value_i=2**i
        select_bits_i =(index&select_value_i)//select_value_i
        graph[i,i]=1
        for j in range(i+1,bits_n):
            select_value_j=2**j
            select_bits_j =(index&select_value_j)//select_value_j
            select_bits   =select_bits_i&select_bits_j
            if abs(sum(ps*select_bits)-pb[i]*pb[j])>1e-12:
                graph[bits_n-i-1,bits_n-j-1]=1
                graph[bits_n-j-1,bits_n-i-1]=1
    #################################################################    
    #最后看y测量独立性
    g=rx(np.pi/2)
    for i in range(1,bits_n):
        g=G(ry(-np.pi/2),g)
    ps     = get_p(g@space)           #统计图
    pb     = [0]*bits_n               #每位是1态的概率
    index  = np.array(range(len(space)))
    for i in range(bits_n):           #获取每位是1态的概率
        pass
        select_value=2**i
        select_bits =(index&select_value)//select_value
        pb[i]       =sum(ps*select_bits)
    for i in range(bits_n):
        select_value_i=2**i
        select_bits_i =(index&select_value_i)//select_value_i
        graph[i,i]=1
        for j in range(i+1,bits_n):
            select_value_j=2**j
            select_bits_j =(index&select_value_j)//select_value_j
            select_bits   =select_bits_i&select_bits_j
            if abs(sum(ps*select_bits)-pb[i]*pb[j])>1e-12:
                graph[bits_n-i-1,bits_n-j-1]=1
                graph[bits_n-j-1,bits_n-i-1]=1
    #################################################################
    cnt=0
    for i in range(bits_n):
        for j in range(i,bits_n):
            if visit[j]==0 and graph[i,j]==1:
                visit[j]=1
                order[cnt]=j
                cnt=cnt+1
    return order
#根据无遮挡序列重新排序各个位并返回重新排序后的态空间
def reorder(space,order):
    sorder=list(range(len(order)))
    for i in range(len(order)):
        if sorder[i]!=order[i]:
            j=sorder.index(order[i])
            space=swap(i,j,space)
            sorder[i],sorder[j]=sorder[j],sorder[i]
    return space

演示代码

隐形传态协议仿真

#隐形传态实验
qubits=zero_qubits(3)#初始全为零态
[qubits[0]]=rand_qubits(1)#随机设置待发送bit
space=get_space(qubits)#获取态空间
print('原始状态')
for v in get_v(space):
    print(v)
#创建纠缠资源
space=GN(1,3,H)@space
space=C1(1,2,3,X)@space
print('创建纠缠资源')
for v in get_v(space):
    print(v)
#发送端
space=C1(0,1,3,X)@space
space=GN(0,3,H)@space
print('发送端测量前的状态')
for v in get_v(space):
    print(v)
print('无遮挡序列：',get_unmusk_order(space))
#接收端
space=C1(1,2,3,X)@space
space=C1(0,2,3,Z)@space
print('接收端恢复后的状态')
for v in get_v(space):
    print(v)

输出：
原始状态
[[ 0.218+0.j ]
[-0.035+0.975j]]
[[1.+0.j]
[0.+0.j]]
[[1.+0.j]
[0.+0.j]]
创建纠缠资源
[[ 0.218+0.j ]
[-0.035+0.975j]]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
发送端测量前的状态
[[ 0.109+0.j ]
[-0.017+0.488j]
[-0.017+0.488j]
[ 0.109+0.j ]
[ 0.109+0.j ]
[ 0.017-0.488j]
[ 0.017-0.488j]
[ 0.109+0.j ]]
无遮挡序列： [0, 1, 2]
接收端恢复后的状态
[[0.707+0.j]
[0.707+0.j]]
[[0.707+0.j]
[0.707+0.j]]
[[ 0.218+0.j ]
[-0.035+0.975j]]

去遮挡演示

#去遮挡演示
#获取全零的态空间
qubits=zero_qubits(4)
space=get_space(qubits)
#使第0位与第3位纠缠
space=GN(0,4,H)@space
space=C1(0,3,4,X)@space
#因为“遮挡”导致分解失败
print('去遮挡前')
for v in get_v(space):
    print(v)
#去除遮挡
order=get_unmusk_order(space)
print('无遮挡序列：',order)
print('去遮挡后')
for v in get_v(reorder(space,order)):
    print(v)

输出：
去遮挡前
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]
[0. +0.j]]
无遮挡序列： [0, 3, 1, 2]
去遮挡后
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
[[1.+0.j]
[0.+0.j]]
[[1.+0.j]
[0.+0.j]]

性能测试

#性能测试
qubits=zero_qubits(10)
space=get_space(qubits)
for i in range(5):
    space=GN(i,10,H)@space
    space=C1(i,i+5,10,X)@space
order=get_unmusk_order(space)
print('无遮挡序列：',order)
for v in get_v(reorder(space,order)):
    print(v)

耗时约1s输出：
无遮挡序列： [0, 5, 1, 6, 2, 7, 3, 8, 4, 9]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]
[[0.707+0.j]
[0. +0.j]
[0. +0.j]
[0.707+0.j]]

其它说明

如果安装了qiskit可以用如下方式在bloch球上画出本文的量子态

import numpy as np
from qiskit.visualization import plot_bloch_vector
#交互
# %matplotlib notebook
#嵌入
%matplotlib inline

def ab2xyz(state):
    #去除全局相位
    if abs(state[0,0])>0:
        state=state*(abs(state[0,0])/state[0,0])
    #转化为正常的坐标(x:横，y:复，z:纵)
    x = np.dot(np.array([[0,1  ]]),state)[0,0].real
    y = np.dot(np.array([[0,-1j]]),state)[0,0].real
    z = np.dot(np.array([[1,0  ]]),state)[0,0].real
    #用二倍角公式转化成bloch sphere坐标
    x = 2*x*z
    y = 2*y*z
    z = 2*z*z-1
    return x,y,z

演示代码：

final_state = np.array([[1],[0]])
plot_bloch_vector(ab2xyz(final_state))

输出：
请添加图片描述
演示代码：

final_state = np.array([[1],[0]])
final_state = rx(np.pi/2)@final_state
plot_bloch_vector(ab2xyz(final_state))

输出：
请添加图片描述

Nervous_46216553

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录