Skip to content

Commit

Permalink
Merge branch 'main' into mpo
Browse files Browse the repository at this point in the history
  • Loading branch information
bartandrews authored Nov 13, 2024
2 parents b9cf3a1 + 87c1e76 commit f6c7fe6
Show file tree
Hide file tree
Showing 3 changed files with 140 additions and 131 deletions.
19 changes: 4 additions & 15 deletions python/ffsim/linalg/double_factorized_decomposition.py
Original file line number Diff line number Diff line change
Expand Up @@ -532,11 +532,7 @@ def double_factorized_t2(
nocc, _, nvrt, _ = t2_amplitudes.shape
norb = nocc + nvrt

occ, vrt = np.meshgrid(range(nocc), range(nvrt), indexing="ij")
occ = occ.reshape(-1)
vrt = vrt.reshape(-1)
t2_mat = t2_amplitudes[occ[:, None], occ[None, :], vrt[:, None], vrt[None, :]]

t2_mat = t2_amplitudes.transpose(0, 2, 1, 3).reshape(nocc * nvrt, nocc * nvrt)
outer_eigs, outer_vecs = _truncated_eigh(t2_mat, tol=tol, max_vecs=max_vecs)
n_vecs = len(outer_eigs)

Expand Down Expand Up @@ -643,16 +639,9 @@ def reconstruct_t2_alpha_beta(
nocc_a, nocc_b, nvrt_a, nvrt_b = t2_amplitudes.shape
norb = nocc_a + nvrt_a

occ_a, vrt_a = np.meshgrid(range(nocc_a), range(nvrt_a), indexing="ij")
occ_b, vrt_b = np.meshgrid(range(nocc_b), range(nvrt_b), indexing="ij")
occ_a = occ_a.reshape(-1)
vrt_a = vrt_a.reshape(-1)
occ_b = occ_b.reshape(-1)
vrt_b = vrt_b.reshape(-1)
t2_mat = t2_amplitudes[
occ_a[:, None], occ_b[None, :], vrt_a[:, None], vrt_b[None, :]
]

t2_mat = t2_amplitudes.transpose(0, 2, 1, 3).reshape(
nocc_a * nvrt_a, nocc_b * nvrt_b
)
left_vecs, singular_vals, right_vecs = _truncated_svd(
t2_mat, tol=tol, max_vecs=max_vecs
)
Expand Down
136 changes: 136 additions & 0 deletions tests/python/states/sample_slater_test.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
# (C) Copyright IBM 2024.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Tests for sampling Slater determinants."""

from __future__ import annotations

import itertools

import numpy as np
import pytest

import ffsim
from ffsim.states.bitstring import BitstringType


@pytest.mark.parametrize(
"norb, nelec, bitstring_type",
[
(norb, nelec, bitstring_type)
for (norb, nelec), bitstring_type in itertools.product(
ffsim.testing.generate_norb_nelec(range(1, 5)), BitstringType
)
],
)
def test_sample_slater_determinant_spinful(
norb: int, nelec: tuple[int, int], bitstring_type: BitstringType
):
"""Test sample Slater determinant, spinful."""
rng = np.random.default_rng(1234)
shots = 1000
for _ in range(min(2, ffsim.dim(norb, nelec))):
rotation_a = ffsim.random.random_unitary(norb, seed=rng)
rotation_b = ffsim.random.random_unitary(norb, seed=rng)
occupied_orbitals = ffsim.testing.random_occupied_orbitals(
norb, nelec, seed=rng
)
rdm_a, rdm_b = ffsim.slater_determinant_rdms(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
vec = ffsim.slater_determinant(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b),
norb,
nelec,
shots=shots,
bitstring_type=bitstring_type,
seed=rng,
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


@pytest.mark.parametrize(
"norb, nelec, bitstring_type",
[
(norb, nelec, bitstring_type)
for (norb, nelec), bitstring_type in itertools.product(
ffsim.testing.generate_norb_nocc(range(1, 5)), BitstringType
)
],
)
def test_sample_slater_determinant_spinless(
norb: int, nelec: int, bitstring_type: BitstringType
):
"""Test sample Slater determinant, spinless."""
rng = np.random.default_rng(1234)
shots = 1000
rotation = ffsim.random.random_unitary(norb, seed=rng)
for occupied_orbitals in itertools.combinations(range(norb), nelec):
rdm = ffsim.slater_determinant_rdms(norb, occupied_orbitals, rotation, rank=1)
vec = ffsim.slater_determinant(norb, occupied_orbitals, rotation)
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
rdm, norb, nelec, shots=shots, bitstring_type=bitstring_type, seed=rng
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


def test_sample_slater_determinant_large():
"""Test sample Slater determinant for a larger number of orbitals."""
norb = 6
nelec = (3, 2)

rng = np.random.default_rng(1234)
shots = 5000
rotation_a = ffsim.random.random_unitary(norb, seed=rng)
rotation_b = ffsim.random.random_unitary(norb, seed=rng)
occupied_orbitals = ((0, 2, 3), (2, 4))
rdm_a, rdm_b = ffsim.slater_determinant_rdms(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
vec = ffsim.slater_determinant(norb, occupied_orbitals, (rotation_a, rotation_b))
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b), norb, nelec, shots=shots, seed=rng
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


def test_sample_slater_determinant_restrict():
"""Test sample Slater determinant with subset of orbitals."""
norb = 8
nelec = (4, 3)

shots = 10
occupied_orbitals = ((0, 2, 3, 5), (2, 3, 4))
rdm_a, rdm_b = ffsim.slater_determinant_rdms(norb, occupied_orbitals)
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b), norb, nelec, orbs=([1, 2, 5], [3, 4, 5]), shots=shots
)
assert samples == ["011110"] * 10
116 changes: 0 additions & 116 deletions tests/python/states/slater_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
import pytest

import ffsim
from ffsim.states.bitstring import BitstringType


@pytest.mark.parametrize("norb, nelec", ffsim.testing.generate_norb_nocc([4, 5]))
Expand Down Expand Up @@ -69,118 +68,3 @@ def test_slater_determinant_amplitudes_spinful(norb: int, nelec: tuple[int, int]
orbital_rotation=(orb_rot_a, orb_rot_b),
)
ffsim.testing.assert_allclose_up_to_global_phase(actual, expected)


@pytest.mark.parametrize(
"norb, nelec, bitstring_type",
[
(norb, nelec, bitstring_type)
for (norb, nelec), bitstring_type in itertools.product(
ffsim.testing.generate_norb_nelec(range(1, 5)), BitstringType
)
],
)
def test_sample_slater_determinant_spinful(
norb: int, nelec: tuple[int, int], bitstring_type: BitstringType
):
"""Test sample Slater determinant, spinful."""
rng = np.random.default_rng(1234)
shots = 1000
for _ in range(min(2, ffsim.dim(norb, nelec))):
rotation_a = ffsim.random.random_unitary(norb, seed=rng)
rotation_b = ffsim.random.random_unitary(norb, seed=rng)
occupied_orbitals = ffsim.testing.random_occupied_orbitals(
norb, nelec, seed=rng
)
rdm_a, rdm_b = ffsim.slater_determinant_rdms(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
vec = ffsim.slater_determinant(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b),
norb,
nelec,
shots=shots,
bitstring_type=bitstring_type,
seed=rng,
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


@pytest.mark.parametrize(
"norb, nelec, bitstring_type",
[
(norb, nelec, bitstring_type)
for (norb, nelec), bitstring_type in itertools.product(
ffsim.testing.generate_norb_nocc(range(1, 5)), BitstringType
)
],
)
def test_sample_slater_determinant_spinless(
norb: int, nelec: int, bitstring_type: BitstringType
):
"""Test sample Slater determinant, spinless."""
rng = np.random.default_rng(1234)
shots = 1000
rotation = ffsim.random.random_unitary(norb, seed=rng)
for occupied_orbitals in itertools.combinations(range(norb), nelec):
rdm = ffsim.slater_determinant_rdms(norb, occupied_orbitals, rotation, rank=1)
vec = ffsim.slater_determinant(norb, occupied_orbitals, rotation)
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
rdm, norb, nelec, shots=shots, bitstring_type=bitstring_type, seed=rng
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


def test_sample_slater_determinant_large():
"""Test sample Slater determinant for a larger number of orbitals."""
norb = 6
nelec = (3, 2)

rng = np.random.default_rng(1234)
shots = 5000
rotation_a = ffsim.random.random_unitary(norb, seed=rng)
rotation_b = ffsim.random.random_unitary(norb, seed=rng)
occupied_orbitals = ((0, 2, 3), (2, 4))
rdm_a, rdm_b = ffsim.slater_determinant_rdms(
norb, occupied_orbitals, (rotation_a, rotation_b)
)
vec = ffsim.slater_determinant(norb, occupied_orbitals, (rotation_a, rotation_b))
test_distribution = np.abs(vec) ** 2
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b), norb, nelec, shots=shots, seed=rng
)
addresses = ffsim.strings_to_addresses(samples, norb, nelec)
indices, counts = np.unique(addresses, return_counts=True)
assert np.sum(counts) == shots
empirical_distribution = np.zeros(ffsim.dim(norb, nelec), dtype=float)
empirical_distribution[indices] = counts / shots
assert np.sum(np.sqrt(test_distribution * empirical_distribution)) > 0.99


def test_sample_slater_determinant_restrict():
"""Test sample Slater determinant with subset of orbitals."""
norb = 8
nelec = (4, 3)

shots = 10
occupied_orbitals = ((0, 2, 3, 5), (2, 3, 4))
rdm_a, rdm_b = ffsim.slater_determinant_rdms(norb, occupied_orbitals)
samples = ffsim.sample_slater_determinant(
(rdm_a, rdm_b), norb, nelec, orbs=([1, 2, 5], [3, 4, 5]), shots=shots
)
assert samples == ["011110"] * 10

0 comments on commit f6c7fe6

Please sign in to comment.