Test that connections that should be are actually permutations

[inducer]

This adds a test that currently fails: Restriction to the boundary and opposite-face swap fall into the "interpolation matrix" code path here rather than the permutation matrix one, even though the matrices should be permutations. There's some crappy thresholding right before that to decide which is which. This is an efficiency problem, as operations that should merely be indirect memory access actually end up consuming considerable flops.

cc @lukeolson

[alexfikl]

I played around a bit with this locally and it can pass the tests if the tolerance in here gets pushed up a bit.

Most of the fails seem to be from the opposite face connection, which also seems to have its own tolerance that's quite a bit larger. Matching them up a bit seems to make most of the tests pass.

[inducer]

Thanks for taking a look! Maybe be it'll be enough to twiddle tolerances to get this to go. I agree that the opposite-face-match tolerance and the permutation-matrix one should be roughly matched up, although they're really for different things. (nodes vs interpolation matrix entries) I'm a bit hesitant to dial down the Gauss-Newton tolerance in opposite-face matching without looking at the history of why it got dialed up that much in the first place. :)

[alexfikl]

I'm a bit hesitant to dial down the Gauss-Newton tolerance in opposite-face matching without looking at the history of why it got dialed up that much in the first place. :)

I tried to dial it down to 1.0e+3 and it started failing pretty fast at higher orders. Even 1.0e+4 started failing for orders 7-8. To be precise, there's an assert in there that started yelling, not the actual Newton iteration (scratch that, the convergence check failed too pretty much at the same time)

[inducer]

Is it promising to try and build the interpolation matrix in quad precision? (float128) It's a one-time cost...

(Possible catch: not all CPU architectures give us float128 with numpy, notably Power9, aka one of the main CEESD execution platforms 😢 )

[alexfikl]

Is it promising to try and build the interpolation matrix in quad precision? (float128) It's a one-time cost...

I don't know if that will work, since numpy.linalg doesn't like float128 from what I remember, so modepy.resampling_matrix will raise. Do you mean just casting it after the fact to do the row sums?

EDIT: Just checked and

>>> x = np.random.rand(5, 5).astype(np.float128)
>>> la.solve(x, np.random.rand(5))
Traceback (most recent call last):
TypeError: array type float128 is unsupported in linalg

[inducer]

Thanks for checking! AFAIR, all we need is linalg.solve. This crappy row echelon factorizer might get us halfway there, we'd only need some triangular solves. (We know the matrices are invertible, so RREF=LU).

(One of the tri solves could be copied from here: https://relate.cs.illinois.edu/course/cs450-s19/f/demos/upload/linear_systems/Coding%20back-substitution.html)

[inducer]

I just spent a bit more time here, and I learned a few things.

• First, when trying to find the permutation between unit nodes, there's kind of no point to constructing the interpolation matrix only to match it to the identity. We might as well compare the nodes directly. It's O(n^2) instead of O(n^3) to boot. That seems to help some cases. (193fabf)
• After that modification, the only cases still failing the assertion that the restriction to the boundary is a permutation are the higher-order warp-and-blend nodes. And the test might be correct: 2D facial restrictions of the 3D warp-and-blend nodes don't seem to be equivalent to 2D warp-and-blend nodes. (They are for all the other node types.) As a result, I think we need to introduce a separate element group type for 2D restrictions of 3D warp-and-blend.
• Plenty of cases still fail the assertion that opposite-face-swap is a permutation. I picked 3D second-order equispaced as (what I thought was) the most egregious example to fail this test. What I found was that the Gauss-Newton objective function seems to have multiple minima, one at the "wrong answer", and one at the right answer. I added a debug snippet to the opposite-face finding (hence the flake8 failures) that proves the existence of those two minima:
  

  meshmode/meshmode/discretization/connection/opposite_face.py

  Lines 225 to 264 in a3c689b

  	if 0:
  	near_zero = np.abs(src_unit_nodes) < 1e-5
  	near_one = np.abs(src_unit_nodes-1) < 1e-5
  	near_mone = np.abs(src_unit_nodes+1) < 1e-5
  	reasonable = near_zero | near_one | near_mone
  	if np.any(~reasonable):
  	unreasonable_els, unreasonable_nodes = np.where(
  	np.any(~reasonable, axis=0))
  	mapped = apply_map(src_unit_nodes)
  	# print(src_unit_nodes[~reasonable])
  	urel = unreasonable_els[0]
  	tgt_urel = tgt_bdry_nodes[:, urel]
  	src_urel = src_bdry_nodes[:, urel]
  	dist_vecs = tgt_urel.reshape(3, -1, 1) - src_urel.reshape(3, 1, -1)
  	dists = np.sqrt(np.sum(dist_vecs**2, axis=0))
  	mat = (dists < 1e-12)
  	from_indices = np.array([np.where(row)[0][0] for row in mat])
  	sunodes = src_bdry_discr.groups[i_src_grp].unit_nodes
  	print("GOOD UNIT NODES")
  	print(sunodes[:, from_indices])
  	print("BAD UNIT NODES")
  	print(src_unit_nodes[:, urel])
  	print()
  	fixed_src_unit_nodes = src_unit_nodes.copy()
  	fixed_src_unit_nodes[:, urel] = sunodes[:, from_indices]
  	mapped_fixed = apply_map(fixed_src_unit_nodes)
  	print("FIXED RESIDUAL")
  	print(np.max(np.abs(mapped_fixed[:, urel]
  	- tgt_bdry_nodes[:, urel])))
  	print("UNFIXED RESIDUAL")
  	print(np.max(np.abs(mapped[:, urel] - tgt_bdry_nodes[:, urel])))
  	pu.db

  
  To see this, you need to change the if 0 to if 1 and run

  pycl test_connection.py 'test_bdry_restriction_is_permutation(_acf, PolynomialEquidistantSimplexGroupFactory, 3, 2)'
  

  which will print something like

  GOOD UNIT NODES
  [[-1. -1. -1.  0.  0.  1.]
   [-1.  0.  1. -1.  0. -1.]]
  BAD UNIT NODES
  [[-1.00000000e+00 -1.00000000e+00 -9.84191183e-01  1.65458682e-16
    -1.11094593e-15  1.00000000e+00]
   [-1.00000000e+00  5.28164206e-16  1.00000000e+00 -1.00000000e+00
    -2.07264860e-15 -1.00000000e+00]]
  
  FIXED RESIDUAL
  1.1102230246251565e-15
  UNFIXED RESIDUAL
  1.3322676295501878e-15
  

  Note the -0.984 hanging out in a set of unit nodes that should only consist of -1s, 0s, and 1s. (Remember, equispaced, order 2!)

I'm out of steam for the night, but I hope this will help get us closer to tracking this down.

(cc @nchristensen since we had just discussed this in our meeting)

[nchristensen]

Interesting. Thanks for looking!

[inducer]

#225 is a version of this that's supposed to be mergeable right away, but with some tests still xfail'd.

[inducer]

#225 and #226 should resolve the practical impact of this. There remains the mystery of why Gauss-Newton had two minima; I'd really like to know why. My plan is to update this after those two merge to leave this as an experimentation ground to investigate that question.

[inducer]

Continued in #228.