Region diff fix

[krooken]

The purpose of the function region_diff() is to remove a set of polytopes (a region, the subtrahend) from one polytope (the minuend). The result may be a single polytope if the result is convex, or multiple polytopes (a single region).

The function region_diff() does not produce the correct result for certain input. For instance, when the minuend is a unit square and the subtrahend is two squares with side 0.5 positioned at [(0, 0), (0.5, 0.5)] and [(0.5, 0.5), (1, 1)], the result produced by region_diff() incorrectly includes the square [(0, 0), (0.5, 0.5)]. The issue is the code that prunes the search tree when none of the remaining polytopes removes anything from the current set of hyperplanes.

The algorithm is a depth first search algorithm implemented in a while loop. It has a set of hyperplanes which is initialized to the minuend. It then works by branching on the hyperplanes of the subtrahend. At each branchin point the set of hyperplanes is partitioned into two by one of the hyperplanes of the subtrahend.

A recursive version of the algorithm is described in Algorithm 2 in:

Mato Baotic "Polytopic computations in constrained optimal control"
Automatika, 2009

When none of the remaining polytopes intersects the current set of hyperplanes, some things hold:

• the current set of hyperplanes is a non-empty partition and is obviously part of the difference
• the algorithm is done in the current branch
• the other partition is empty, since the intersection between the current set of hyperplanes and the subtrahend is empty

This PR:

• adds tests to assert correctness of some simple cases
• adds one test that fails without this fix
• tries to improve the readability a bit
• adds comments
• fixes the failing test case

[murrayrm]

It would be great to get a review of this from @tichakornw or @johnyf if they have time, since they would be able to spot errors better than I could.

[krooken]

The more the merrier. I thought it was quite difficult to find a solution, so it might be a good idea to do some iterations before accepting this PR.

[johnyf]

Thank you for the changes, and for organizing them as separate commits.

In summary, except for:

• the algorithm change in commit e7036dc, which I currently do not understand (please see the relevant review comments for more details), but which does make a failing test pass (this newly introduced test, namely test_checker_region_diff in script tests/polytope_test.py, appears to me correct, which indicates that there is an error in the existing implementation of the function polytope.polytope.region_diff), and
• a new call introduced in commit 66da390 to function polytope.polytope.reduce which appears to increase worst-case complexity, be motivated as an optimization, and for which it is unclear whether the average-case complexity is improved, or what the performance effects are, and
• the code comments introduced in commit 765a7a9, some of which I need to study more,

the other changes look good to me, subject to the modifications needed that are
described in my review comments (and other suggestions in the review comments).

About the change in commit e7036dc and the related test that demonstrates an error before this change (function test_checker_region_diff in script tests/polytope_test.py, which would better be introduced in a separate commit just before commit e7036dc, please see the review comments at function test_checker_region_diff for more details), it may be better to separate them into another pull request.

This would also align with a comment above (#67 (comment)), which I understand to refer to the change of commit e7036dc.

About commit messages, I would like to recommend the following:

• to use categorizing prefixes for commit titles ("TST: ", "REF: ", "BUG: ", etc.)

• as already true for some commits in this pull request, to use verbs in the imperative mood in the commit title, as recommended for the Linux kernel:

  Describe your changes in imperative mood, e.g. "make xyzzy do frotz" instead of "[This patch] makes xyzzy do frotz" or "[I] changed xyzzy to do frotz", as if you are giving orders to the codebase to change its behaviour.

  This practice allows more to fit in the title's bounded width, forms shorter commit messages, and more readable commit titles (see also this).

• to use Markdown in commit messages

• to characterize items in the commit body, for example as "function", "class", etc. This is really the same recommendation as for naming expressions in mathematical writing, for example "the variable x" instead of "x", or "the matrix A", instead of "A". Relevant book: "Mathematical Writing" by Donald E. Knuth, Tracy Larrabee, and Paul M. Roberts, AMS, 1996.
  For example, possible modifications to the body of the message of commit cfba0be are:

  • "the builtin function sum" or "the expression sum(counter)"
  • "an instance of numpy.ndarray"
    (in this rephrasing, the package numpy is also mentioned to identify class ndarray)
  • "the function numpy.sum"

• except for the prefix, to use lowercase for the commit title (exceptions may be names of persons, or analogous cases)

For example, applying the above to commit titles in this pull request could result in the following titles:

• "TST: test function polytope.polytope.region_diff"
• "MAI: call function reduce within region_diff"
  ("MAI" is used here because this change is not a refactoring.)
• "REF: change to pythonic indexing"
• "REF: remove unused return values"
• "REF: improve readability of function region_diff"
• "REF: remove superfluous if statement"
• "BUG: correct indexing error"
  ("BUG" signifies that this is a correction)
• "BUG: correct algorithm of function region_diff"
• "DOC: comment in region_diff"

The prefx "REF:" explicitly communicates the claim that the changes in the commit result in equivalent code, i.e., that they leave the implemented algorithm unaffected.

More information about commit prefixes (list of the above and other prefixes, their meaning, and origin in numpy development practices), as well as other resources, including links to recommendations about writing commit messages, can be found in tulip's documentation: https://github.com/tulip-control/tulip-control/blob/bb004422b575dccea8d19c33acfeb04b37c62a5a/doc/dev_guide.rst#advice

[johnyf]

This reference could be expanded: Automatika, Vol. 50, No. 3--4, pp. 119--134, 2009.

I do not know whether there exists a DOI URL to the paper, but a URL to a PDF may be relevant. I found a PDF of this journal paper at: https://hrcak.srce.hr/file/71840, but not at the publisher: https://www.tandfonline.com/loi/taut20, because the list of issues available online at the publisher apparently starts at Vol. 51 (2010), and this paper is in Vol. 50 (2009). The journal Print ISSN: 0005-1144 on the publisher's website matches the ISSN noted within the paper's PDF. I found the journal's contents from earlier years at: https://hrcak.srce.hr, and the issues containing this paper at: https://hrcak.srce.hr/index.php?show=toc&id_broj=3787. The paper title is also listed on what seems to be the author's website. There is also a Scopus entry.

[johnyf]

The second name of the author is Baotić. Earlier in the same file the name is written in UTF8:

polytope/polytope/polytope.py

Line 41 in 0ccd0f6

# M. Kvasnica, P. Grieder and M. Baotić,

so it can be written also here using UTF8 (alternatively it could be written using LaTeX as \'{c}).

[johnyf]

About blank lines within functions (indeed there already are such blank lines throughout this test script), PEP 8 recommends:

Use blank lines in functions, sparingly, to indicate logical sections.

I would recommend not inserting blank lines within functions. I agree that there may be areas within functions where some vertical separation could be helpful to the reader, after all, PEP 20 recommends:

Sparse is better than dense.

but I think that vertical separation would better be obtained by:

• inserting comments, in particular "sectioning" comments, i.e., replacing the blank line with a (preferably) one-line comment that fits as a title to describe the block of lines of code that follow.

  Since this block may already contain comment lines, one way to emphasize the "title" comment is to add a preceding line with an empty comment. What the empty line achieves is to introduce vertical separation, while at the same time serving as a visual guide of the indentation depth within the function. Maintaining a continuous line of indentation is visually important to both connect the different parts of a function, and to form one uninterrupted block of nonempty lines for each function. So this approach would increase readability, following PEP 20:

  Readability counts.

• inserting blank comments (same reasons as for blank comments in the previous item).

An example from the same test script:

polytope/tests/polytope_test.py

Lines 398 to 418 in 0ccd0f6

	def test_lpsolve():
	# Ensure same API for both `scipy` and `cvxopt`.
	# Ensured by the different testing configurations.
	# Could change `polytope.polytope.default_solver` to
	# achieve the same result, when `cvxopt.glpk` is present.
	#
	# 2-D example
	c = np.array([1, 1], dtype=float)
	A = np.array([[-1, 0], [0, -1]], dtype=float)
	b = np.array([1, 1], dtype=float)
	res = solvers.lpsolve(c, A, b)
	x = res['x']
	assert x.ndim == 1, x.ndim
	assert x.shape == (2,), x.shape
	#
	# 1-D example
	c, A, b = example_1d()
	res = solvers.lpsolve(c, A, b)
	x = res['x']
	assert x.ndim == 1, x.ndim
	assert x.shape == (1,), x.shape

This comment applies throughout the code changes.

[johnyf]

I am in favor of keeping the first decimal digit, for example (0.0, 1.0). This also matches the string representation of float in Python str(1.0) == '1.0'.

I think that the result is more readable, and also according to PEP 20:

Explicit is better than implicit.

This comment applies throughout the code changes.

[johnyf]

Please use keyword-value syntax for calling the method str.format, for example:

assert diff == remaining, '{d}, {r}'.format(
    d=diff, r=remaining)

This makes the code more readable, and less fragile to changes.

Also, I strongly recommend printing the values in the same order they appear in the assertion's expression: first diff, then remaining. This is more readable. Swapping them can easily confuse the reader of the program's output.

I am also in favor of just using tuples in this case (of an assertion that tests equality):

assert diff == remaining, (diff, remaining)

This is succinct, since the source code line that contains the assert statement will be printed for the user to read when the assertion fails.

Finally, I am in favor of single quotes as string delimiters: 'foo' instead of "foo". Single quotes introduce less visual clutter than double quotes (PEP 20: "Sparse is better than dense."), so the code becomes more readable (PEP 20: "Readability counts."). Also, this style for string delimiters is already used in the rest of this test script.

This comment applies throughout the code changes.

[johnyf]

It is not clear to me why this change (i.e., commit e7036dc) is justified (other than the change of the test result for the function test_checker_region_diff in the script tests/polytope_test.py from failing in commit cfba0be to passing in commit e7036dc, specifically in the last assertion within that function: https://github.com/tulip-control/polytope/pull/67/files#diff-67aa94aa544367d400036bb409ff79a2d399229092060d7ac8d28b8f1c38499bR559) (I should note that the newly introduced test test_checker_region_diff appears correct to me, which indicates that the existing implementation of function polytope.polytope.region_diff contains an error). This is a big change that appears to change the algorithm.

For example, in the removed lines (lines 2079 and 2081--2096 on the left) there are a break statement and a return statement. In contrast, the replacement code remains within the function's body.

Having said this, I will continue trying to understand how this change affects the algorithm.

[krooken]

I will try to address the concerns in several parts. Let's first look at the return statement.

For the program to jump to the else statement on the removed line 2091 the value of counter[level] must necessarily be strictly greater than 1 (and usually more). This is so because mi can only have positive elements. There is only one line on which counter[level] can be increased past 1, and that is on removed line 2104:

counter[level] = counter[level] + 1

Now, it is not sufficient for counter[level] to be strictly greater than 1 for the program to reach the else statement on line 2091. counter[level] must also be strictly greater than mi[level]. This can never happen at this point in the program. Right after the increment of counter[level] on line 2104 counter[level] is reset to 0 if counter[level] is greater than mi[level]:

counter[level] = counter[level] + 1
if counter[level] <= mi[level]:
    ...
else:
    counter[level] = 0

Hence, the else statement on line 2091 cannot be reached and the return statement is dead code. I think this is why the indexing error on line 2093 has gone unnoticed.

[krooken]

Let's also look at the for loop on removed line 2082. The iteration variable jj is never used, so we effectively have the following code on removed lines 2079-2096:

level = level - 1
res = union(res, Polytope(A[INDICES, :], B[INDICES]), False)
INDICES[-1] -= M
INDICES = np.hstack([
    INDICES,
    beg_mi[level] + counter[level] + M
])

[krooken]

When running the test test_checker_region_diff in file polytope_test.py the values of some variables when the while loop on line 2071 is entered:

INDICES = [0, 1, 2, 3]
mi = [2, 2]

The content of INDICES represents the four sides of the minuend polytope. The values of mi tells us that the first and the second (ordered by decreasing Cheby ball radius) polytope in the subtrahend region has two hyperplanes each that are not found in the minuend polytope.

The halfplanes are indexed like this:
0-3: The four halfplanes of the minuend
4, 5: The two halfplanes of the first subtrahend
6, 7: The two halfplanes of the second subtrahend
8, 9: The complements of each of the halfplanes of the first subtrahend
10, 11: The complements of each of the halfplanes of the second subtrahend

First pass through while loop:

level == 0
counter == [0, 0]

The program is at first pass at level 0. The two hyperplanes (indices 4 and 5) of the largest subtrahend (happens to be the lower left box) is intersected with the minuend. The intersection is non-empty, so there is a need to see what is left of the minuend. So the program branches on the first hyperplane (index 4) of the first subtrahend, and it looks at the side of the hyperplane on which the first subtrahend is not (index 8). The first halfplane happens to be x < 0.5, so now the program looks at x > 0.5 to see whether this part of the minuend is non-empty. It is, so the level is increased.

Second pass through while loop:

level == 1
counter == [1, 0]
INDICES = [0, 1, 2, 3, 8]

The program is at first pass at level 1. The two hyperplanes (indices 6 and 7) of the remaining subtrahend (the upper right box) is intersected with the right part of the minuend (right of 0.5). The intersection is non-empty, so there is a need to see what is left of the right part of the minuend. The program branches on the first hyperplane (index 6) of the remaining subtrahend, and it looks on the side on which the subtrahend is not (index 10). This halfplane happens to be x > 0.5, so now the program looks at x < 0.5 to see whether this part of the right part of the minuend is non-empty. It is not, so the level remains at the current value.

Third pass through while loop:

level == 1
counter == [1, 1]
INDICES == [0, 1, 2, 3, 8, 10]

The program is at second pass at level 1. Now it needs to look at the other side of the first hyperplane of the second subtrahend (x > 0.5, index 6), and on the side of the second hyperplane (index 7) on which the second subtrahend is not (y < 0.5, index 11), to see if this part of the minuend is non-empty. It is non-empty, so the program can continue. But since the program now is at the highest possible level (there are no more subtrahends to subtract) the current set of halfplanes (indices [1, 2, 3, 8, 6, 11], representing the polytope 0.5 < x < 1.0 and ´0.0 < y < 0.5`) is a subset of the resulting region difference. Therefore, it is added to the result.

Fourth pass through while loop:

level == 1
counter == [1, 2]
INDICES == [0, 1, 2, 3, 8, 6, 11]

The program is done at this level, so it pops out of this level by decrementing the level and resetting the counter and the indices to a previous state, and then advances the state to the next. The next state should check x < 0.5 and y > 0.5 (indices 4 and 9). This part of the minuend is non-empty, so the program increments the level to see whether there are any subtrahends left to subtract.

Fifth pass through the while loop:

level == 1
counter == [2, 0]
INDICES == [0, 1, 2, 3, 4, 9]

The program is at the first pass at level 1, but the current set of halfplanes does not intersect with the remaining subtrahend. Therefore the program enters the if statement on removed line 2078. Here things start to go wrong. The conclusion at this point, that the current set of indices represents a subset of the resulting difference, is correct. But then the program should back out of the current level (nothing more here to do) and advance the state. However, what happens is that the level is decremented and the program will next check 0.0 < x < 0.5 and 0.0 < y < 0.5 (indices 0, 1, 2, 3, 4, 5, and 10). This is a non-empty set, and it does not intersect with the second subtrahend, so the first subtrhend is added to the result. This is of course not good.

The problem is that the program does not pop out of the current level in a correct way. The correct way to pop a level is implemented in the else statement on line 2112. To get there, INDICES needs to represent an empty set, and counter[level] must indicate that the program is done at the current level (counter[level] must be at least mi[level]).

[johnyf]

I agree with the reasoning for this change. However, this change modifies code statements, but it appears in commit 765a7a9, which is titled "Adding comments to region_diff".

I strongly recommend splitting commit 765a7a9 into a commit with this change to code statements, and another commit with the changes to code comments.

Also, about the change itself, I would recommend:

if len(mi) > 1:
    csum = np.cumsum(np.hstack([0, mi[:-1]]))

[johnyf]

This line is too long (109 characters). Please shorten this line. According to PEP 8:

Limit all lines to a maximum of 79 characters.
...
The limits are chosen to avoid wrapping in editors with the window width set to 80, even if the tool places a marker glyph in the final column when wrapping lines.

(I actually write up to 80 characters, but usually I make lines shorter than that.)

This comment applies to several other comment lines in commit 765a7a9.

[johnyf]

This line is too long (110 characters). Please shorten this line. According to PEP 8:

Limit all lines to a maximum of 79 characters.
...
The limits are chosen to avoid wrapping in editors with the window width set to 80, even if the tool places a marker glyph in the final column when wrapping lines.

(I actually write up to 80 characters, but usually I make lines shorter than that.)

This comment applies to one more comment line in commit e7036dc.

[johnyf]

This assertion fails before commit e7036dc, and passes with that commit. So it may be worth splitting the code of function test_checker_region_diff into a separate commit, and shifting that commit (which would introduce only the function test_checker_region_diff) to become the commit preceding commit e7036dc, i.e., to group these commits along the history, because they are relevant.

This rewriting of history can be obtained with two calls to git rebase: the first git rebase edits commit 76773b3, by amending it to omit this function, and committing the omitted lines as a subsequent commit, then git rebase --continue to complete the first rebase, and the second git rebase shifts the newly created commit to just before (the commit on the rewritten history that corresponds to) commit e7036dc.

I lean towards keeping commits that add tests that demonstrate errors separate from the commits that correct the errors, and placing the former before the latter in git history, so that the failure can be confirmed by running the tests with the repository at the former commit.

Still, it may be better to both introduce tests that demonstrates the error and correct the error within a single commit. I do that too. It may be more helpful as documentation for changes that are difficult to understand without an example that demonstrates the error.

[krooken]

Nice! Thank you for the review comments. I will start looking at them this week, but it will take some time to go through them.

I agree that the change of the algorithm in commit e7036dc is quite big and that it is difficult to understand how the change affects the search tree. If it was implemented with recursion it would be simpler to understand and to correct, but such a solution would potentially use a lot of stack memory.

I went through the while loop several times when I was debugging. I will try to write down step-by-step how the old implementation fails and how the new version fixes it.

[johnyf]

Indeed, I agree that a recursive implementation would be more readable, and that in CPython it would also be more limited, due to the potential for raising a RecursionError.

[johnyf]

Yet another observation in support of using the imperative mood when titling commit messages: the changes to the LaTeX kernel and to LaTeX packages in the LaTeX news have such titles, for example: https://www.latex-project.org/news/latex2e-news/ltnews30.pdf. This is an addendum to my comment above (#67 (review)).

[krooken]

I think all review comments should be addressed now. I tried to explain how and why the algorithm failed before and how the fix solves it, but the code change is still the same as before.

[johnyf]

Thank you for revising the changes, I will review the revised changes.

[johnyf]

A note that while reviewing the changes of this pull request, I had encountered issues in a function in the call graph of the function region_diff (I think it is the function polytope.polytope.reduce). The issues may require re-implementing the function reduce, and possibly also other functions in the module polytope.polytope.

These issues caused errors in the tests of tulip. A re-implementation of parts of the module polytope.polytope may be needed, to avoid these errors. The current code (i.e., in polytope == 0.2.3) in module polytope.polytope can compute results that are in error, without raising an exception.