Двухрежимная функция Вигнера в питоне - PullRequest
Новая
       19

Двухрежимная функция Вигнера в питоне

0 голосов
/ 31 октября 2019

Я попытался настроить функцию Вигнера, чтобы она обрабатывала два состояния режима, особенно для итеративного метода.

Однако размер массива, который выдает мой вывод, слишком велик, и я не знаю, почему? Именно тогда, когда я пытаюсь вычислить логарифмическую отрицательность Вигнера, используя ее, интегралы получаются в виде массивов, а не единичных значений.

Ниже приведен код и описание того, что предполагается сделать: 

`import numpy as np
from scipy import (zeros, array, arange, exp, real, conj, pi,
               copy, sqrt, meshgrid, size, polyval, fliplr, conjugate,
               cos, sin)
import scipy.sparse as sp
import scipy.fftpack as ft
import scipy.linalg as la
from scipy.special import genlaguerre
from scipy.special import binom
from scipy.special import sph_harm

from qutip.qobj import Qobj, isket, isoper
from qutip.states import ket2dm
from qutip.parallel import parfor
from qutip.utilities import clebsch
from scipy.special import factorial
from qutip.cy.sparse_utils import _csr_get_diag
from qutip import *

def wigner2(psi, xvec1, yvec1, xvec2, yvec2, method='iterative', g=np.sqrt(2)):
    """Wigner function for a state vector or density matrix at points
     `xvec1 + i * yvec1` `xvec2 + i * yvec2`

Параметры

state: qobj Вектор состояния или матрица плотности.

xvec1: массив_подобных x-координат, по которым вычисляется функция Вигнера.

yvec1: массив_подобныхy-координаты, по которым вычисляется функция Вигнера.

xvec2: массив_подобных x-координат для вычисления функции Вигнера.

yvec2: массив_подобных y-координат для вычисления функции Вигнера.

g: float Коэффициент масштабирования для a = 0.5 * g * (x + iy), по умолчанию g = sqrt(2).

method: string {'iterative'} Выберите метод «iterative», где «iterative» использует итерационный метод дляоценить функции Вигнера для матриц плотности: математика: |m><n|. Метод «итеративный» используется по умолчанию и в целом рекомендуется, но метод «Лагерра» более эффективен для матриц очень разреженной плотности (например, суперпозиции состояний Фока в большом гильбертовом пространстве). Метод 'fft' является предпочтительным методом для работы с матрицами плотности, которые имеют большое количество возбуждений (> ~ 50).

Возвращает

W: массив Значения, представляющие функцию Вигнера, рассчитанную поуказанный диапазон [xvec1, yvec1] и [xvec2, yvec2] "" "

    if not (psi.type == 'ket' or psi.type == 'oper' or psi.type == 'bra'):
        raise TypeError('Input state is not a valid operator.')
    if psi.type == 'ket' or psi.type == 'bra':
        rho = ket2dm(psi)   #always use  density matrix
    else:
        rho = psi
    if method == 'iterative':
        return _wigner2_iterative(rho, xvec1, yvec1, xvec2, yvec2, g)
    else:
        raise TypeError(
            "method must be 'iterative'")


def _wigner2_iterative(rho, xvec1, yvec1, xvec2, yvec2, g=np.sqrt(2)):

" "" Использование итерационного метода для оценки функций Вигнера для состояния Фока: математика: |mp><nq| ВигнерФункция рассчитывается следующим образом: математика: W = \sum_{mpnq} \\rho_{mpnq} W_{mpnq}, где: математика: W_{mpnq} - функция Вигнера для матрицы плотности: математика: |mp><nq|. В этой реализации для каждой строки m * p Wlist содержит функции Вигнера Wlist = [0, ..., W_mpmp, ..., W_mpnq]. Как только вычисляется одна функция Вигнера W_mpnq, соответствующий вклад добавляется к общей функции Вигнера, взвешенной соответствующим элементом в матрице плотности: math: rho_{mpnq}. "" "

    M1 = np.prod(ptrace(rho, 0).shape[0])
    M2 = np.prod(ptrace(rho, 1).shape[0])
    M = np.prod(rho.shape[0])  
    X1, Y1, X2, Y2 = np.meshgrid(xvec1, yvec1, xvec2, yvec2)
    A1 = 0.5 * g * (X1 + 1.0j * Y1 + 0 * X2 + 0 * Y2)
    A2 = 0.5 * g * (0 * X1 + 0 * Y1 + X2 + 1.0j * Y2)


    Wlist1 = array([zeros(np.shape(A1), dtype=complex) for k in range(M)])

    Wlist2 = array([zeros(np.shape(A2), dtype=complex) for k in range(M)])


    W = real(rho[0, 0]) * real(Wlist1[0] * Wlist2[0])


    for m in range(0,M1):
        if m==0:
            Wlist1[0] = exp(-2.0 * abs(A1) ** 2) / (pi) 
        else:
            Wlist1[m] = ((2.0 * A1 * Wlist1[m - 1]) / sqrt(m))
        for n in range(0, M2):
            if n==0:
                Wlist2[0] = exp(-2.0 * abs(A2) ** 2) / (pi) 
            else:
                Wlist2[n] = ((2.0 * A2 * Wlist2[n - 1]) / sqrt(n))
            if m != 0 and n != 0:
                W += 2 * real(rho[0, m * M2 + n] * Wlist1[m] * Wlist2[n])

    for p in range(0, M1):
        temp1 = copy(Wlist1[m])
        temp2 = copy(Wlist2[n])
        if p==0:
            Wlist1[p] = exp(-2.0 * abs(A1) ** 2) / (pi) 
        else:
            Wlist1[p] = ((2.0 * conj(A1) * temp1 -sqrt(p) * Wlist1[p-1]) / sqrt(p))
            for q in range(0, M2):
                if q==0:
                    Wlist2[q] = exp(-2.0 * abs(A2) ** 2) / (pi) 
            else:
                Wlist2[q] = ((2.0 * conj(A2) * temp2 - sqrt(q) * Wlist1[q - 1]) / sqrt(q))                
                W += 2 * real(rho[p * M2 + q, p * M2 + q] * Wlist1[p] * Wlist2[q])

            if p != 0 and q !=0:
                for k in range(p + 1, M1):
                    temp3 = (2 * A1 * Wlist1[k-1] - sqrt(k) * temp1) / sqrt(k)
                    temp1 = copy(Wlist1[k])
                    Wlist1[k] = temp3
                    for l in range(q +1, M2):
                        temp4 = (2 * A2 * Wlist2[l-1] - sqrt(l) * temp2) / sqrt(l)
                        temp2 = copy(Wlist2[l])
                        Wlist2[l] = temp4
                        W += 2 * real(rho[p * M2 + q, k *M2 +l] * Wlist1[k] * Wlist2[l])

    return 0.5 * W * g **2'
...