Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Windows fix for numba types #389

Merged
merged 14 commits into from
Apr 26, 2024
25 changes: 12 additions & 13 deletions thewalrus/internal_modes/fock_density_matrices.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,13 +79,13 @@ def _density_matrix_single_mode(
x = np.array(x)
Ax = nb_ix(A, x, x) # A[np.ix_(x, x)]

edge_reps = np.array([half_c, half_c, 1] + list(pattern))
edge_reps = np.array([half_c, half_c, 1] + pattern)
n_edges = 3 + K * len(pattern)

assert n_edges == Ax.shape[0] // 2 == 3 + K * (M - 1)

N_max = 2 * edge_reps.sum()
N_fixed = 2 * np.sum(pattern)
N_fixed = 2 * sum(pattern)

steps = np.prod(edge_reps + 1)

Expand Down Expand Up @@ -130,7 +130,9 @@ def _density_matrix_single_mode(
haf_arr += haf_arr_new

rho = (
(-1) ** pattern.sum() * haf_arr / (np.sqrt(np.linalg.det(Q).real) * np.prod(fact[pattern]))
(-1) ** sum(pattern)
* haf_arr
/ (np.sqrt(np.linalg.det(Q).real) * np.prod(fact[np.array(pattern)]))
)

for n in range(cutoff):
Expand Down Expand Up @@ -197,7 +199,7 @@ def density_matrix_single_mode(
K = cov.shape[0] // (2 * M)
if not set(list(pattern.keys())).issubset(set(list(np.arange(M)))):
raise ValueError("Keys of pattern must correspond to all but one spatial mode")
N_nums = np.array(list(pattern.values()))
N_nums = list(pattern.values())
HM = list(set(list(np.arange(M))).difference(list(pattern.keys())))[0]
if LO_overlap is not None:
if not K == LO_overlap.shape[0]:
Expand Down Expand Up @@ -225,19 +227,18 @@ def density_matrix_single_mode(
return vals
if method in ["non-recursive", "diagonals"]:
cov = project_onto_local_oscillator(cov, M, LO_overlap=LO_overlap, hbar=hbar)
num_modes = len(cov) // 2
A = Amat(cov)
Q = Qmat(cov)
pref = 1 / np.sqrt(np.linalg.det(Q).real)
blocks = np.arange(K * M).reshape([M, K])
dm = np.zeros([cutoff, cutoff], dtype=np.complex128)
num_modes = M * K
block_size = K
if method == "non-recursive":
num_modes = M * K
block_size = K
for i in range(cutoff):
for j in range(i + 1):
if (i - j) % 2 == 0:
patt_long = list((j,) + tuple(N_nums) + ((i - j) // 2,))
patt_long = [j] + N_nums + [(i - j) // 2]
new_blocks = np.concatenate((blocks, np.array([K + blocks[-1]])), axis=0)
perm = (
list(range(num_modes))
Expand All @@ -260,11 +261,9 @@ def density_matrix_single_mode(
dm[j, i] = 0
else:
for i in range(cutoff):
patt_long = (i,) + tuple(N_nums)
dm[i, i] = (
pref
* haf_blocked(A, blocks=blocks, repeats=patt_long)
/ np.prod(factorial(patt_long))
patt_long = [i] + N_nums
dm[i, i] = pref * np.real(
haf_blocked(A, blocks=blocks, repeats=patt_long) / np.prod(factorial(patt_long))
nquesada marked this conversation as resolved.
Show resolved Hide resolved
)
if normalize:
dm = dm / np.trace(dm)
Expand Down
Loading