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.

Sign up for GitHub

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

Jump to bottom

Feature/fix repeated open contours #112

Merged
merged 4 commits into from
Nov 17, 2023
Merged

Feature/fix repeated open contours #112

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1d25895
prune repeated contours when open
piskunow Nov 17, 2023
e9a6c8c
fix space
piskunow Nov 17, 2023
38ff882
update tutorial with open and periodic cases
piskunow Nov 17, 2023
166abb3
fix typing
piskunow Nov 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
460 changes: 155 additions & 305 deletions docs/tutorial.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

57 changes: 29 additions & 28 deletions src/fermi_contours/marching_squares.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,8 @@
"""Marching Squares module."""

from typing import Any
import logging
from typing import Callable
from typing import Optional
from typing import Union
from warnings import warn

import numpy as np
import numpy.typing as npt
Expand Down Expand Up @@ -117,8 +115,9 @@ def __call__(self, level: float = 0) -> list[LPFloat]:
contour_paths: list of lists of pairs of floats.
Each list has numerical interpolated points along the path.
"""
_, contour_paths = self._find_contours(level)
return contour_paths
contour_cells, contour_paths = self._find_contours(level)
_, pruned_paths = self._check_repeated(contour_cells, contour_paths)
return pruned_paths

@property
def grid_points(self) -> tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]]:
Expand Down Expand Up @@ -224,9 +223,9 @@ def _find_contours(self, level: float) -> tuple[list[LPInt], list[LPFloat]]:
try:
d_ij = marching_step(cells[ij], self.func, middle_k, d_ij)
except RuntimeError:
warn("Saddle point not resolved.")
logging.debug("Saddle point not resolved.")
if self.func is None:
warn(
logging.debug(
"Saddle point not resolved because 'func' is not provided."
)
break
Expand Down Expand Up @@ -272,7 +271,9 @@ def _find_contours(self, level: float) -> tuple[list[LPInt], list[LPFloat]]:
contour_paths.append(single_path)
break
else:
warn(f"Stepping outside the initial path with cell {ij}.")
logging.debug(
"Stepping outside the initial path with cell %s", ij
)

if (0 <= next_i < n_x) and (0 <= next_j < n_y):
pass
Expand All @@ -292,26 +293,26 @@ def _find_contours(self, level: float) -> tuple[list[LPInt], list[LPFloat]]:

def _check_repeated(
self, contours_cells: list[LPInt], contour_paths: list[LPFloat]
) -> list[LPFloat]:
# check for repeated cells and return only the largest
digested_contour_cells: list[Any] = []
digested_contour_paths: list[Any] = []
for c_cells, v_cells in zip(contours_cells, contour_paths):
# check that this does not belong to digest
updated = False
for digest, paths in zip(
digested_contour_cells,
digested_contour_paths,
):
if not set(c_cells).isdisjoint(digest):
digest.update(c_cells)
updated = True
paths.update(v_cells)
if not updated:
digested_contour_cells.append(set(c_cells))
digested_contour_paths.append(v_cells)

return [list(d.values()) for d in digested_contour_paths]
) -> tuple[list[LPInt], list[LPFloat]]:
# open contours may start several times, prune them to keep the largest
# for each path
pruned_cell_list: list[LPInt] = []
pruned_path_list: list[LPFloat] = []
for contour, path in zip(contours_cells, contour_paths):
# replace subsets in pruned list
for idx, pruned_path in enumerate(pruned_path_list):
if set(path).issuperset(pruned_path):
pruned_cell_list[idx] = contour
pruned_path_list[idx] = path
# add missing sets from the pruned list
is_subset = any(
set(path).issubset(pruned_path) for pruned_path in pruned_path_list
)
if not is_subset:
pruned_cell_list.append(contour)
pruned_path_list.append(path)

return pruned_cell_list, pruned_path_list


def marching_cell_values(
Expand Down
98 changes: 97 additions & 1 deletion tests/test_ms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,82 @@ def surface(x: float, y: float) -> float:
(8, 5),
(8, 4),
]
val_contour_paths_dict = zip(val_contours_cells, val_contour_paths)


level_open = 1.1

val_contour_paths_open = [
[
(-0.5263157894736843, 0.9004328254847646),
(-0.6989414121724066, 0.7777777777777777),
(-0.736842105263158, 0.741871858007592),
(-0.8843425113710203, 0.5555555555555554),
(-0.9473684210526316, 0.4361380424746075),
(-0.9907748538011696, 0.33333333333333326),
(-1.0376884340480832, 0.11111111111111116),
(-1.0376884340480832, -0.11111111111111116),
(-0.9907748538011696, -0.33333333333333337),
(-0.9473684210526316, -0.4361380424746075),
(-0.8843425113710202, -0.5555555555555556),
(-0.736842105263158, -0.741871858007592),
(-0.6989414121724064, -0.7777777777777778),
(-0.5263157894736843, -0.9004328254847646),
(-0.31611842105263155, -1.0),
],
[
(0.5263157894736841, -0.9004328254847647),
(0.6989414121724065, -0.7777777777777778),
(0.7368421052631575, -0.7418718580075925),
(0.8843425113710202, -0.5555555555555556),
(0.9473684210526314, -0.4361380424746081),
(0.9907748538011696, -0.33333333333333337),
(1.0376884340480832, -0.11111111111111116),
(1.0376884340480832, 0.11111111111111116),
(0.9907748538011696, 0.33333333333333326),
(0.9473684210526314, 0.436138042474608),
(0.8843425113710204, 0.5555555555555554),
(0.7368421052631575, 0.7418718580075926),
(0.6989414121724066, 0.7777777777777777),
(0.5263157894736841, 0.9004328254847647),
(0.3161184210526316, 1.0),
],
]
val_contours_cells_open = [
[
(7, 8),
(6, 8),
(6, 7),
(5, 7),
(5, 6),
(4, 6),
(4, 5),
(4, 4),
(4, 3),
(4, 2),
(5, 2),
(5, 1),
(6, 1),
(6, 0),
(7, 0),
],
[
(11, 0),
(12, 0),
(12, 1),
(13, 1),
(13, 2),
(14, 2),
(14, 3),
(14, 4),
(14, 5),
(14, 6),
(13, 6),
(13, 7),
(12, 7),
(12, 8),
(11, 8),
],
]


def test_single_contour() -> None:
Expand All @@ -54,3 +129,24 @@ def test_single_contour() -> None:
assert_allclose(contours_cells[0], val_contours_cells)
assert_allclose(contour_paths[0], contour)
assert_allclose(val_contour_paths, contour)


def test_open_contours() -> None:
"""Test closed single contour."""
squares = MarchingSquares(
func=surface,
res=(20, 10),
bounds=((-2, 2), (-1, 1)),
)

contours_cells, contour_paths = squares._find_contours(level=level_open)
pruned_cells, pruned_paths = squares._check_repeated(contours_cells, contour_paths)

contours = squares(level=level_open)

assert_allclose(pruned_cells[0], val_contours_cells_open[0])
assert_allclose(pruned_cells[1], val_contours_cells_open[1])
assert_allclose(pruned_paths[0], contours[0])
assert_allclose(pruned_paths[1], contours[1])
assert_allclose(val_contour_paths_open[0], contours[0])
assert_allclose(val_contour_paths_open[1], contours[1])
Loading