QISKIT 機能色々(TIPS)

主に自分用のメモです。
コマンドの列挙ですので、同時に動かない機能があります。

• 状態ベクトル系では、回路にmeasureがあるとバグる。（正確には、射影測定後の状態ベクトルが得られる)
• 古典レジスタが存在すると使えない機能がある

install

pip install qiskit pylatexenc

import

#必要なモジュールのインポート
from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute, Aer
from qiskit.qasm import pi
from qiskit.tools.visualization import plot_histogram, circuit_drawer, plot_gate_map, plot_circuit_layout
from qiskit import IBMQ, transpile
from qiskit.providers.ibmq.managed import IBMQJobManager
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.quantum_info import Statevector
from qiskit.aqua.utils import tensorproduct
from qiskit.quantum_info.analysis import average_data
from qiskit.aqua import QuantumInstance
from qiskit.aqua.operators import PauliExpectation, CircuitSampler, StateFn
from qiskit.aqua.operators import X, Y, Z, I
from qiskit.circuit.library import RYGate
from qiskit.aqua.operators.state_fns import CircuitStateFn
from qiskit.circuit.random import random_circuit
from qiskit.circuit.add_control import add_control
from qiskit.providers.aer import QasmSimulator

import numpy as np

IBMQロード

## "MY_API_TOKEN" に、上で取得したtokenを入れる
from qiskit import IBMQ
IBMQ.save_account('XXXXXXXX',overwrite = True)

# least busy な バックエンドを知る
from qiskit.providers.ibmq import least_busy

# 自分のアカウント情報をloadする。（あらかじめ IBMQ.save_account を実行しておく必要がある. 複数のアカウントを使い分ける時はここで行う)
provider = IBMQ.load_account()

# 自分のアカウントで使用できるバックエンドを見る
provider.backends()

# Run our circuit on the least busy backend. Monitor the execution of the job in the queue
from qiskit.tools.monitor import job_monitor

オブザーバブルの作り方

# Observableを作る　(方法その1)
pauli0 = Pauli(label='II')
pauli1 = Pauli(label='IZ')
obs_op = Operator(pauli0) + Operator(pauli1) #O = II+IZ
obs1 = obs_op.data # from operator to matrix
print(obs1)
# Observableを作る　(方法その2)
I_mat = np.array([[1, 0], [0, 1]])
Z_mat = np.array([[1, 0], [0, -1]])
obs2 = tensorproduct(I_mat,I_mat) + tensorproduct(I_mat,Z_mat) # matrix

# Observableを作る　(方法その3)
obs3 = np.array([[2,0,0,0],[0,0,0,0],[0,0,2,0],[0,0,0,0]])

# Observableを作る　(方法その4)
obs4 =  1*I^I + 1*I^Z

# Observableを作る　(方法その5)
 from qiskit.aqua.operators import WeightedPauliOperator
 pauli_dict = {
      'paulis': [{"coeff": {"imag": 0.0, "real": 1}, "label": "II"},
                {"coeff": {"imag": 0.0, "real": 1}, "label": "IZ"},                             
                ]
  }
obs5 = WeightedPauliOperator.from_dict(pauli_dict).to_opflow()

ちなみに、方法4は若干の注意が必要です。
obs4 = 2*I^I + 2*I^Z はOKですが、
obs4 = 2*(I^I + I^Z) はだめです。
(エラー無く通りますが、期待値計算の際に定数倍が無視されます）

なんだかんだ、汎用性が高いのは方法5だと思います。
obs5.to_matrix() で配列形式にも出来ます。

回路準備

#　回路の定義
nqubits = 2 #qubits
ncbits = 1  #classical bits
qr = QuantumRegister(nqubits)
cr0 = ClassicalRegister(ncbits)
cr1 = ClassicalRegister(ncbits)
circuit = QuantumCircuit(qr,cr0,cr1) 
# レジスタに興味ない場合は QuantumCircuit(nqubits)や QuantumCircuit(nqubits,ncbits) でOK
#circuit.draw('mpl')

回路記述いろいろ

# 回路記述

# Define an arbitrary unitary gate
U_mat = np.array([[0,1],[1,0]])

#Note : this operation is invalied in ibmq real quantum devices
#circuit.reset(0) #reset 0th qubit as |0>
#circuit.reset(1) #reset 1st qubit as |0>
#circuit.initialize([0, 1],0) #set 0th qubit as [0,1] = |1>
circuit.h(0)
circuit.h(1)
#circuit.x(0).c_if(cr1,1)　# Note : Only 1 classical bit register can be used in c_if.
#circuit.unitary(U_mat,0)
#circuit.x([0,1]) #add X-gate from 0th to 1st qubit.
#circuit.cx(0,1) #Controlled-NOT. The controlled and target bit is 0th and 1st one, respectively.
#circuit.rz(2,0) ## RZ(a)=exp(-i*a/2*Z)に注意. rz(rotation_angle, qbit_label)
#circuit.hamiltonian(obs1,10,[0,1]) #Hamitonian evaluation ?
psi = CircuitStateFn(circuit) # get final state vector. Note:The classical register shouldn't be used.
#circuit.measure_all() #Note : don't add measurement when you use statevector_simulator
circuit.draw('mpl')

psi = CircuitStateFn(circuit) を使う場合は、
-古典レジスタを使わない
-measureを使わない
が必要になります。

任意制御ゲートの作成

theta_01 = np.arctan(1)
Ugate_01 = RYGate(2*theta_01)
Ugate_01_ctrl = add_control(Ugate_01, num_ctrl_qubits=2, label='CU_01', ctrl_state='01')

def eigen_rot(circ):
    circ.append(Ugate_01_ctrl,[0,1,3])

回路の途中の状態ベクトルをメモる

回路の途中の状態ベクトルを取り出すことが出来ます。
この方法はシミュレーター限定ですが、statevector simulatorに限らず、qasm simulatorでも使えます。
回路の任意の位置で　snapshot という操作を入れます。（ゲートのような感覚です）
circuit.snapshot('my_sv',snapshot_type='statevector') # get a snapshot of present statevector
snapshotはいくつでも入れられます。
次に、回路を実行します。
results = execute(circuit, backend, shots=1).result()
結果の変数resultsの中身に、以下の呪文でアクセスします。
statevec = results.data()['snapshots']['statevector']['my_sv'][k] # 2^nqubits complex vector |ψ>
ここで、kはsnapshotを複数とった場合に何番目を取り出すかを意味します。
出力結果は複素ベクトル（状態ベクトル）になっていますので、期待値計算は簡単です。
expectation = np.real(np.vdot(statevec,np.dot(obs_mat,statevec))) # Real[<ψ|O|ψ>]
obs_matはオブザーバブルの行列表現です。

c_ifの使い方

2021/01/02追記
c_ifは1bit古典レジスタしか引数に取れません。以前は複数bitが出来たようです。
複数bitでも実行可能に見えますが、実際には何も操作されなかったりします。ひどい。。
例えば、

は正常に動きます("11")が、

は正常に動きません("10")。なんでや！　たぶん最初の古典ビットだけがc_ifの対象ぽいですね。

また、複数古典ビットを一気に渡す
cr = ClassicalRegister(2), circuit.x(0).c_if(cr,3) も、無効です。以前は出来た模様。

1 classical bit register をたくさん定義するのは面倒です。
circuit.x(0).c_if(cr[0],1) と書きたい。その場合は以下のようにします。

qr = QuantumRegister(2)
circuit = QuantumCircuit(qr) 
cr = [ ClassicalRegister(1) for _ in range(2) ]
for register in cr:
    circuit.add_register( register )
circuit.x([1])
circuit.measure(1,1)
circuit.x(0).c_if(cr[1],1)
circuit.measure(0,0)

また、c_ifは実機では動かないそうです。以下の5-1を見てください。
https://qiskit.org/textbook/ch-algorithms/teleportation.html#testing

回路を調べる

# depth of circuit
circuit.depth()

## decompose into hardware gates
circuit.decompose().draw('mpl')

backendの定義

# backendの定義
#backend = Aer.get_backend('qasm_simulator')
#backend = Aer.get_backend('statevector_simulator')
#backend = Aer.get_backend('unitary_simulator')
#backend = provider.get_backend('ibmq_santiago')

transpile関係

#ibmq information

#transpile
circuit_transpiled = transpile(circuit, backend=provider.get_backend('ibmq_santiago'), optimization_level=3)
plot_circuit_layout(circuit_transpiled, backend=provider.get_backend('ibmq_santiago'))

IBM実機プロパティを取る

# T1/T2 info.
backend = provider.get_backend('ibmq_santiago')
props = backend.properties()
def describe_qubit(qubit, properties):
    """Print a string describing some of reported properties of the given qubit."""

    # Conversion factors from standard SI units
    us = 1e6
    ns = 1e9
    GHz = 1e-9

    print("Qubit {0} has a \n"
          "  - T1 time of {1} microseconds\n"
          "  - T2 time of {2} microseconds\n".format(
              qubit,
              properties.t1(qubit) * us,
              properties.t2(qubit) * us))
              #properties.gate_error('u2', qubit),
              #properties.gate_length('u2', qubit) * ns,
              #properties.frequency(qubit) * GHz))

describe_qubit(0, props)
describe_qubit(1, props)

qubitアサイン手動設定

# Virtual -> physical
#    0    ->    3
#    1    ->    4
#    2    ->    2

circuit_opt = transpile(circuit, backend=provider.get_backend('ibmqx2'), initial_layout=[3, 4, 2])
plot_circuit_layout(circuit_opt, backend=provider.get_backend('ibmqx2'))

#execute関連

IBMQ認証

"MY_API_TOKEN" に、IBMQ HPで取得したtokenを入れる

from qiskit import IBMQ
IBMQ.save_account('MY_API_TOKEN',overwrite = True)

# 自分のアカウント情報をloadする。（あらかじめ IBMQ.save_account を実行しておく必要がある. 複数のアカウントを使い分ける時はここで行う)
provider = IBMQ.load_account()

# 自分のアカウントで使用できるバックエンドを見る
print(provider.backends())

# 実行（実行状況のモニタアリ）
from qiskit.tools.monitor import job_monitor
job= execute(circuit, backend=provider.get_backend('ibmq_qasm_simulator'), shots=8192) #ibmq
results = job.result()
counts=results.get_counts()
plot_histogram(counts)

単一jobのexecute

# execute on qasm simulator
results = execute(circuit, backend=Aer.get_backend('qasm_simulator'), shots=1024).result()

# execute on ibmq
#job= execute(circuit, backend=provider.get_backend('ibmq_santiago'), shots=1024)
#job_monitor(job, interval=2)
#results = job.result()
#counts = results.get_counts()

# execute on state vector simulator
#results = execute(circuit, backend=Aer.get_backend('statevector_simulator'),shots=1).result()
#results.get_statevector()

過去のjobの回収

指定したマシンで直近10個の実行結果のIDを取る

[x.job_id() for x in provider.get_backend('ibmq_qasm_simulator').jobs(limit=10)]

指定したID"hogehoge"の結果を回収する

provider.get_backend('ibmq_qasm_simulator').retrieve_job("hogehoge")

jobの実行上限を確認

backend_job_limit = backend.job_limit()
remain_jobs = backend.remaining_jobs_count()
maximum_jobs = backend_job_limit.maximum_jobs
print('Job is {:} / {:}'.format(remain_jobs, maximum_jobs))

Job is 5 / 5

error mitigation

from qiskit.aqua import QuantumInstance
# Import measurement calibration functions
from qiskit.ignis.mitigation.measurement import (complete_meas_cal, tensored_meas_cal,CompleteMeasFitter, TensoredMeasFitter)
qi_noise_model_ibmq = QuantumInstance(backend=provider.get_backend('ibmqx2'), shots=8000, measurement_error_mitigation_cls=CompleteMeasFitter, measurement_error_mitigation_shots=8000)
job=qi_noise_model_ibmq.execute(circuit_opt)
results = job # q_instance.execute の戻り値は、既に.resultなので注意。
counts=results.get_counts()
plot_histogram(counts)

noise model を使ったexecute

#Noise model
from qiskit.providers.aer.noise import NoiseModel
backend_dev = provider.get_backend('ibmqx2')
coupling_map = backend_dev.configuration().coupling_map
noise_model = NoiseModel.from_backend(backend_dev.properties())
q_instance = QuantumInstance(backend=Aer.get_backend("qasm_simulator"), shots=8192, noise_model=noise_model,  coupling_map=coupling_map, measurement_error_mitigation_cls=CompleteMeasFitter, cals_matrix_refresh_period=30)
job=q_instance.execute(circuit)
results = job
counts=results.get_counts()
plot_histogram(counts)

または、

  backend_dev = provider.get_backend('ibmqx2')
  coupling_map = backend_dev.configuration().coupling_map
  basis_gates =  backend_dev.configuration().basis_gates
  noise_model = NoiseModel.from_backend(backend_dev.properties())
  print(noise_model)  
  results = execute(circs,backend=backend,   noise_model=noise_model, coupling_map=coupling_map, basis_gates=basis_gates).result()  

または、

ibmqx2_simulator = QasmSimulator.from_backend(provider.get_backend('ibmqx2'))
results = execute(circs,backend=ibmqx2_simulator).result()

複数jobの同時execute

 # Build ten circuits.
circs = []
for _ in range(10):
    circs.append(random_circuit(n_qubits=5, depth=4))

results = execute(circs,backend=backend).result()
results.get_counts(5)  # Counts for experiment 5.

または、

 # Build ten circuits.
circs = []
for _ in range(10):
    circs.append(random_circuit(n_qubits=5, depth=4))

job_manager = IBMQJobManager()
job_set_foo = job_manager.run(circs, backend=backend, name='foo')
results = job_set_foo.results()
results.get_counts(5)  # Counts for experiment 5.

IBMQ job manager を使う場合、transpileがすっ飛ばされる。
そのため、特定のゲートの回路しか通らない。
これを回避するには、手動transpileしてから投げればOK。(backend指定transpileは、やや長めの時間がかかります）
circs = transpile(circs, backend)
https://qiskit.org/documentation/stubs/qiskit.providers.ibmq.managed.IBMQJobManager.html

測定結果の加工

results.get_memory(circuit) #全shotの結果を列挙
counts = results.get_counts()
counts.most_frequent() # 頻度最大のbit列を取得
counts.get('00') # '00'の観測された回数を取得

期待値計算

#Calc. expectation

# qasm_simulator or ibmq

#方法その1 , Note : obs should be diagonalized in Z basis
#counts = results.get_counts()
#expectation = average_data(counts, obs)

#方法その2
q_instance = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), shots=8192)
# define the state to sample
measurable_expression = StateFn(obs4,is_measurement=True).compose(psi) 
# convert to expectation value
expectation = PauliExpectation().convert(measurable_expression)  
sampler = CircuitSampler(q_instance).convert(expectation)    # get state sampler (you can also pass the backend directly)
Expectation = sampler.eval().real

# state vector simulator
#result_statevec = Statevector(results.get_statevector())
#expectation = np.real(result_statevec.expectation_value(Operator(obs1)))

plot関連

ヒストグラムplot

plot_histogram(counts)

確率0のcountsも表示する

# zero-padding
for i in range (2**nqubits):
    counts.setdefault(format(i, '0'+str(nqubits)+'b'),0)

plot_histogram(counts, figsize = [20,5])

statevector simulator を使ったときのヒストグラムplot

n = nqubits
results = execute(circuit, backend=Aer.get_backend('statevector_simulator'), shots=1).result()
final_state = results.get_statevector()

## z方向に観測した時の確率分布を求める. (状態ベクトルの各成分の絶対値の二乗=観測確率)
probs = np.abs(final_state**2)
print(probs)
# プロットする
import matplotlib.pyplot as plt
%matplotlib inline

## z方向に射影測定した時に得られる可能性があるビット列
z_basis = [format(i,"b").zfill(n) for i in range(probs.size)]

plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()

Statevectorクラスのメソッドを使っても出来るかもしれません。

回路プロットの横幅変更

改行したくないときはfoldの値を増やす。

circuit.draw('mpl', fold = 30)

更に、 scale = 0.5 等で画像サイズを変えられる。

matplotlibによるグラフ

"""複数のグラフを重ねて描画するプログラム"""
import matplotlib.pyplot as plt

fig, ax = plt.subplots()


c1,c2 = "blue","green"     # 各プロットの色
l1,l2 = "ibmq_santiago","simulator"   # 各ラベル

ax.set_xlabel('time')  # x軸ラベル
ax.set_ylabel('magnetization')  # y軸ラベル
ax.set_title('Evolution of magnetization') # グラフタイトル
# ax.set_aspect('equal') # スケールを揃える
ax.grid()            # 罫線
#ax.set_xlim([-10, 10]) # x方向の描画範囲を指定
#ax.set_ylim([0, 1])    # y方向の描画範囲を指定
ax.plot(x, y, color=c1, label=l1, marker="o")
ax.plot(x_sim, y_sim, color=c2, label=l2)

ax.legend(loc=0)    # 凡例
fig.tight_layout()  # レイアウトの設定
# plt.savefig('hoge.png') # 画像の保存
plt.show()

qiskit.aqua

QAOA

# Weighted matrix. 
w = np.array([
    [0, 1, 0, 1],
    [1, 0, 1, 0],
    [0, 1, 0, 1],
    [1, 0, 1, 0]
])

p = 1
backend = Aer.get_backend('statevector_simulator')
optimizer = COBYLA()
qubitOp, offset = max_cut.get_operator(w)

qaoa = QAOA(qubitOp, optimizer, p)
quantum_instance = QuantumInstance(backend)
result = qaoa.run(quantum_instance)


Statevec = Statevector(result.eigenstate)
counts = Statevec.probabilities_dict()
plot_histogram(counts)

その他

circuit.remove_final_measurement() の意義

動作がよくわかりません。
古典レジスタも一緒に消してしまって、回路として使えなくなってしまうような。

状態ベクトルシミュレーターのときはmeasureを無効化するオプションがほしい

QPEのライブラリ