密度演算子(DensOpクラス)
========================

## 密度演算子の操作

量子状態には純粋状態と混合状態があります。混合状態は、複数の純粋状態が
ある確率でもって出現するアンサンブル{p(i), |phi(i)>} (i=0,1,2,...)とし
て定義されます。このような混合状態も含めて、量子状態を統一的に記述する
便利な記法として「密度演算子」が導入されました。その定義は、密度演算子
をrhoとすると、「rho = ∑_i p(i) |phi(i)><phi(i)|」です。

qlazyには、通常の量子状態の集合から密度演算子をつくり、各種計算する機
能があります。

### 生成・初期化　

DensOpクラスが密度演算子を表します。使用例を以下に示します。

    >>> from qlazy import QState,DensOp
	>>> 
    >>> qs1 = QState(2).h(0).cx(0,1)
    >>> qs2 = QState(2).x(0).z(1)
	>>> 
    >>> de = DensOp(qstate=[qs1,qs2], prob=[0.3,0.7])

２つの量子状態qs1,qs2があったとき、DensOpの引数として、量子状態のリス
トと各々の出現確率を表すリストを与えると、密度演算子のインスタンスが返っ
てきます。量子状態は、何個指定しても良いです（メモリが許す限り）。ただ
し、指定する量子状態の量子ビットはすべて一致していなければなりません。
また、probを省略した場合、各量子状態が等しい確率で混合しているものとし
て密度演算子が計算されます。

numpyの行列（２次元配列）を指定して生成することもできます。使い方は、
以下です。

    >>> import numpy as np
	>>> mat = np.array([[1,0],[0,0]])
    >>> de = DensOp(matrix=mat)
	
行列の次元は2のべき乗である必要があります。

### ゲート演算

密度演算子として表現された量子状態に対してゲート演算することができます。
適用するゲート演算子をU、密度演算子をrhoとしたとき、「U * rho * U^」
という行列計算を実行しています。仕様はQStateクラスの場合と全く同じです。
以下に使用例を示します。

    >>> de.h(0).cx(0,1)
	>>> de.crx(0,1, phase=0.1)
    >>> ...

#### カスタム・ゲートの追加

DensOpクラスを継承することで、自分専用の量子ゲートを簡単に作成・追加す
ることができます。ベル状態を作成する回路をbellメソッドとして、DensOpを
継承したMyDensOpクラスに追加する例を示します。

    >>> class MyDensOp(DensOp):
    >>>     def bell(self, q0, q1):
    >>>         self.h(q0).cx(q0,q1)
    >>>         return self
	>>> 
    >>> de = MyDensOp(qubit_num=2)
	>>> de.bell(0,1)
    >>> ...

これは非常に簡単な例なのであまりご利益を感じないかもしれませんが、大規
模な量子回路を作成したい場合など、便利に使える場面は多いと思います。


### パウリ積の演算

パウリ演算子X,Y,Zのテンソル積を定義して密度演算子に演算することができま
す。パウリ積を扱うために、まず、

    >>> from qlazy import DensOp, PauliProduct
	
のようにPauliProductクラスをimportする必要があります。例えば、3量子ビッ
トの密度演算子deに対して、X2 Y0 Z1というパウリ積を演算したい場合、

	>>> pp = PauliProduct(pauli_str="XYZ", qid=[2,0,1])
	>>> de.operate_pp(pp=pp)
	
のようにoperate_ppメソッドのppオプションにPauliProductのインスタンスを指
定します。制御化されたパウリ積はoperate_ppメソッドのqctrlオプションに制御
量子ビット番号を指定することで実現できます。以下のようにします。

	>>> pp = PauliProduct(pauli_str="XYZ", qid=[2,0,1])
	>>> qs.operate_pp(pp=pp, qctlr=3)

### 量子回路の演算

量子状態に量子回路を演算することができます。量子回路はQCircクラスを使って以下のように用意します。

    >>> from qlazy import QCirc
    >>> qc = QCirc().h(0).cx(0,1)
	>>> qc.show()
    q[0] -H-*-
    q[1] ---X-

これを密度演算子に演算するために'operate_qcirc'メソッドを使います。

    >>> de = DensOp(qubit_num=2)
	>>> de.operate_qcirc(qc)

'operate_pp'と同様に制御ビットを追加することもできます。

	>>> de.operate_qcirc(qc, qctrl=3)

ここで１点注意事項があります。演算できる量子回路はユニタリに限ります。
測定ゲートのような非ユニタリゲートを含むものは演算できません。


### 量子チャネル

密度演算子に対して、いくつかの代表的な量子チャネルを適用することができます。

    >>> de.bit_flip(q, prob=xxx)
    >>> de.phase_flip(q, prob=xxx)
    >>> de.bit_phase_flip(q, prob=xxx)
    >>> de.depolarize(q, prob=xxx)
    >>> de.amp_dump(q, prob=xxx)
    >>> de.phase_dump(q, prob=xxx)

上から順に、「ビット反転」、「位相反転」、「ビット・位相反転」、「分極
解消」、「振幅ダンピング」、「位相ダンピング」です。引数のqは適用した
い量子ビット番号です。probは（雑音）適用の確率を表すパラメータです。各々
の量子チャネルが何を意味するかについては、例えば、ニールセン・チャンの
「量子コンピュータと量子通信３」を参照してください。

### メモリ解放

密度演算子インスタンスのメモリは使用されなくなったら自動的に解放されますが、
明示的に解放したい場合、

    >>> del de

とすれば好きなタイミングで解放することができます。

クラス・メソッド'del_all'を使えば複数の密度演算子のインスタンスを一気に
解放することができます。

	>>> mat = np.array([[0.3, 0.1], [0.1, 0.7]])
    >>> de_0 = DensOp(matrix=mat)
    >>> de_1 = DensOp(matrix=mat).x(0)
    >>> de_2 = DensOp(matrix=mat).h(0)
    >>> ...
    >>> DensOp.del_all(de_0, de_1, de_2)

このとき、引数に指定するのは密度演算子のリストやタプルであっても良いで
すし、それらの入れ子でもOKです。例えば、

    >>> de_A = [de_1, de_2]
    >>> DensOp.del_all(de_0, de_A)
	>>> 
    >>> de_B = [de_3, [de_4, de_5]]
    >>> DensOp.del_all(de_B)

という指定の仕方でも大丈夫です。

### 複製

密度演算子を複製することができます。シミュレータならではの機能です。以下
のようにcloneメソッドを呼び出します。

    >>> de_clone = de.clone()

### リセット

すでに生成済みの密度演算子を破棄することなく、再度初期化して使いたい場
合は、reset()メソッドを使います。

    >>> de.reset()

とすれば、deは強制的に|00...0><00...0|となります。これは、同じ量子回路
を何度も繰り返し適用しながら、何らかの統計量を得たいような場合に便利な
機能です。

また、

    >>> de.reset(qid=[1,5])

のように、量子番号リストを引数で与えると、その番号に対応した量子ビット
のみを強制的に|0>にすることもできます。

## 密度演算子の表示

### 要素の表示

DensOpクラスのshowメソッドを使います。使用例を以下に示します。

    >>> de.show()

例えば、以下のように各行列要素の値が表示されます。

    elm[0][0] = +0.0000+0.0000*i : 0.0000 |
    elm[0][1] = +0.0000+0.0000*i : 0.0000 |
    elm[1][0] = +0.0000+0.0000*i : 0.0000 |
    elm[1][1] = +1.0000+0.0000*i : 1.0000 |+++++++++++
	
部分系の表示も可能です。

    >>> de.show(qid=[q0,q1,...])

例えば、

    >>> de.show(qid=[0,3])

とすると、0番目と3番目の量子ビットに関する密度演算子を表示します。内部
では0番目と3番目の量子ビット以外の量子ビットに関して部分トレースを実行
しています。

### 非ゼロ要素のみ表示

ゼロでない要素だけを表示したい場合、以下のように'nonzero'オプションを使えばできます。

    >>> de.show(nonzero=True)

とすると、

    elm[1][1] = +1.0000+0.0000*i : 1.0000 |+++++++++++

のように非ゼロ要素のみが表示されます。

### 要素の取得

get_elmメソッドで取得できます。

    >>> elm = de.get_elm()
    >>> elm = de.get_elm(id=[1,2])


## 密度演算子に関する計算

### トレース

traceメソッドを使います。使用例を以下に示します。

    >>> tr = de.trace()

密度演算子のトレースは理論的には「常に1である」ということになっている
ので、たいていの場合、あまり意味はないメソッドですが、密度演算子に対し
て、何らかの行列変換を適用しておいて(後述)、そのトレースを計算したい場
合があります。その場合、このメソッドを使うことで、計算できます。

### ２乗のトレース

DensOpクラスのsqtraceメソッドを使います。使用例を以下に示します。

    >>> sqtr = de.sqtrace()

密度演算子を２乗してから、そのトレースをとります。純粋状態の密度演算子
の場合、この値は1になり、混合状態の場合、1以下になることが、理論的にわ
かっています。このメソッドを使う典型的な状況は、純粋状態と混合状態を判
別したいときです。

### 部分トレース

DensOpクラスのpatraceメソッドを使います。使用例を以下に示します。

    >>> de_reduced = de.patrace(qid=[q0,q1,...])

全体系の部分系に対する密度演算子を知りたい場合に「部分トレース」という
処理を施せば良いということがわかっています。その計算を実行します。返却
値は、縮約された(部分系の)密度演算子インスタンスです。

### 部分系

部分トレースをとった残りの部分系を取得します。qidに取得したい部分系に
対応した量子ビット番号リストを指定します。

    >>> de_partial = de.partial(qid=[q0,q1,...])
	
patraceメソッドでは「捨てたい」量子ビット番号を指定するのに対し、
partialメソッドでは「残したい」量子ビット番号を指定します。引数の解釈
の違いを除けば、両者の内部処理は全く同じです。

### テンソル積

２つの密度演算子のテンソル積を計算できます。使用例を以下に示します。

    de_product = de_A.tenspro(de_B)

密度演算子de_Aとde_Bのテンソル積(de_A (X) de_B)を計算し、新たな密度演
算子de_productを生成しています（(X)は丸付きの☓のつもり）。

### 複合状態の生成

compositeメソッドを使います。

    >>> de_com = de.composite(4)

全く同じdeを4つ複合した状態が生成されます。

### 行列の適用

applyメソッドを使います。使用例を以下に示します。

    >>> import numpy as np
    >>> M = np.array([[0.0,1.0],[1.0,0.0]])
    >>> de.apply(matrix=M, qid=[2])  # 2番目の量子ビットにXゲート演算

密度演算子rhoに対して変換行列Mがあったときに、「M * rho * M^」の
計算をします。行列はnumpyの２次元配列として用意しておきます。密度演算
子よりも変換行列のサイズの方が小さくなくてはいけません。また、qidで指定
する引数のリストのサイズをnとしたとき、行列の行数、列数はともに2^nでな
くてはなりません。applyメソッドは、元のインスタンスの内容を変えます。
ご注意ください。

### 和とスカラー積

addメソッド、mulメソッドを使います。

    >>> de.add(de_0)
    >>> de.mul(0.3)

deにde_0を足して、さらに0.3倍する例です。add,mulはもとのインスタンスの
内容を変えます。

### 密度演算子の混合

    >>> de = DensOp.mix(densop=[de1,de2], prob=[0.2,0.8])

量子状態ではなく、複数の密度演算子を用意し、その混合で新たな密度演算子
を生成することができます。

### 確率

    >>> prob = de.probability(povm=[E0,E1], qid=[0,1])
    >>> prob = de.probability(kraus=[M0,M1], qid=[0,1])

kraus演算子(numpy行列)のリスト、またはPOVM演算子(numpy行列)のリストを
引数に指定し、qidで測定対象となる量子ビット番号リストを指定します。各々
の測定演算に対する確率値のリストを返します。qidで指定する引数のリスト
のサイズをnとしたとき、kraus演算子(numpy行列)、またはPOVM演算子の行数、
列数はともに2^nでなくてはなりません。

### 測定（CP-instrument）

    >>> de.instrument(kraus=[M0,M1], qid=[0,1])
    >>> de.instrument(kraus=[M0,M1], qid=[0,1], measured_value=1)

引数krausでKraus演算子(numpy行列)のリスト、引数qidで測定対象となる量子
ビット番号リストを指定します。密度演算子をKraus表現に変更します。これ
は測定を実行するがその結果を忘れた状況に相当します（非選択的）。引数
measured_valueに測定値（Kraus演算子番号）を指定すると、密度演算子が測
定後の状態に変化します（選択的）。ただし、正規化しないので、正規化した
い場合は、以下のようにトレースで確率を求めて、それで割る必要があります。

    >>> prob = de.trace()
    >>> de.mul(factor=1.0/prob)

### フィデリティ（忠実度）

	>>> fid = de1.fidelity(de2)

密度演算子de1とde2のフィデリティを計算します。

### トレース距離

	>>> dis = de1.distance(de2)

密度演算子de1とde2のトレース距離を計算します。

### スペクトル分解

	>>> qstate,prob = de.spectrum()

密度演算子をスペクトル分解します。固有値（確率）のリストと固有ベクトル
（量子状態）のリストをリターンします。ちなみに、

	>>> de1 = DensOp(qstate=qstate, prob=prob)

とやると、スペクトル分解の逆操作となり密度演算子は元に戻ります（計算誤
差を除き、ですが）。

### エントロピー（フォンノイマン・エントロピー）

	>>> ent = de.entropy()

密度演算子deのフォンノイマン・エントロピーを計算します。

### エンタングルメント・エントロピー

	>>> ent_A = de.entropy(qid=[0,1])    # for system A
	>>> ent_B = de.entropy(qid=[2,3,4])  # for system B

全体系に対する部分系のエントロピー（エンタングルメント・エントロピー）
を計算します。entropyメソッドに部分系を表す量子ビット番号リストを引数
として与えます。上の例では、全体系が５量子ビットで記述されている想定の
もとで、それに対する部分系A（量子ビット番号:[0,1]）および部分系B（量子
ビット番号:[2,3,4]）のエンタングルメント・エントロピーを計算しています。
全体系が純粋状態の場合、部分系Aと部分系Bのエンタングルメント・エントロ
ピーの値は一致します。

### 条件付きエントロピー

	>>> ent_cond = de.cond_entropy([0,1],[2,3,4])

部分系B（量子ビット番号:[2,3,4]）が定まったという条件のもとでの、部分
系A（量子ビット番号:[0,1]）の条件付きエントロピーを計算します。
S(A|B)=S(A,B)-S(A)と定義されます。古典情報理論における条件付きエントロ
ピーと違い、負の値になる場合もあります。

### 相互情報量

	>>> mut_info = de.mutual_info([0,1],[2,3,4])

部分系A（量子ビット番号:[0,1]）と部分系B（量子ビット番号:[2,3,4]）の相
互情報量を計算します。I(A:B)=S(A)+S(B)-S(A,B)と定義されます。AとBを入
れ替えても同じ値になります。

### 相対エントロピー

	>>> rel_ent = de1.relative_entropy(de2)
	
密度演算子de1とde2の相対エントロピーを計算します。


以上