Add linear pixel boundaries in silicon sensor

[jamesp-epcc]

This draft pull request adds linear pixel boundary arrays alongside the existing polygonal boundaries. It also makes some changes to the original boundary code to prevent adjacent pixels' boundaries from diverging over time. This involves averaging the distortions read in for shared boundary points at startup, and ensuring that both sides of the outermost boundary are updated when updating pixel distortions. Currently the code keeps both sets of boundaries up-to-date, with the useNewBoundaries flag determining which are used for inside pixel tests.

The whole test suite now passes, including a new test to compare the old and new boundaries. Some of the expected values in the tests had to be modified due to the fixes made to the original code. For a more rigorous test, the checkPixel call in updatePixelDistortions can be uncommented. This will compare the linear and polygonal boundary points every time they change and exit with an error as soon as any discrepancy is found.

Performance tests have not yet been carried out. This will be done once the general approach is approved and the old code is removed.

[rmjarvis]

@jamesp-epcc I'm hoping to get to reviewing this soon. But it's still listed as work in progress and not passing unit tests. Can you get this into shape for review soon? (I.e. passing unit tests, full code coverage)

The failure I can see looks like it's testing something at really too many digits of precision, so can either dial down the accuracy requirement on the test and/or change the value it's comparing against.

[jamesp-epcc]

Apologies, not sure how I missed that CI failure. I have now fixed the initial failure which was due to small differences in the areas due to the changes I've made. I am not sure what's causing the latest failure on macOS as it doesn't look to be obviously related to the code I've changed and I can't reproduce it on Linux. I do have access to a Mac so I will see if I can reproduce it on there.

(I flagged this as Draft because I don't think you'll want to merge it until the old code has been removed, but it would be useful to get feedback at this stage).

[rmjarvis]

The MacOS failure is probably spurious. The GH actions mac platforms aren't super stable, so they have a tendency to seg fault pretty much randomly from time to time. I restarted them. If that still fails, then I'll look into it, but usually that makes the problem go away. (I've seen this in multiple repos, and it's a different place in the tests every time, so AFAICT, it's a bug on their end, not in GalSim.)

[rmjarvis]

Thanks James. I finally finished going through this. I thought it was basically fine until I got to the unit tests. I'm actually pretty concerned about the pixel area of the central pixel changing by as much as it did. I know we discussed this on some telecons a while back, but looking at the actual numbers involved, I'm pretty sure this is introducing an error that is more that the accuracy Craig had verified the old code. So probably not such a great thing to change.

I'm proposing an adjustment to the step where you are currently averaging the values from neighboring pixels. I think instead we should just trust the inner pixel's values and ignore the ones from the outer pixel. Both sides' values were verified by Craig to be accurate, so either one should actually be fine, but the average may not be. And I suspect that the corners especially get the wrong value using this procedure.

There are also various style and commentary suggestions I mentioned as I was going.

[rmjarvis]

Hm. This change is actually a little worrisome. I know we changed how the linear interpolation is done between the old and new, but I'm concerned that new version might not be very accurate. I think Craig verified the old number between the old approximation and a much finer grid with the Poisson solver to significantly better than 0.2% (which is the level of change here for the central pixel).
The maximum value changing by ~8e-5 seems ok. But I'm not so sure about the (0,0) change of 2e-3.

[rmjarvis]

I wonder if the problem is mostly coming at the corners. The corners of pixel (0,0) move basically directly inward. So they are pretty accurate from the (0,0) pixel's own coordinates for that point, as the angle of the 45 degree direction is still pointing exactly at the corner.
For the upper right corner, the (1,1) pixel's lower left direction probably still has a good estimate for this point.
But the (0,1) and (1,0) pixels will not be aiming at the right place for their upper left and lower right directions (respectively), so they will hit different points on the curve that aren't actually at the corner point. Averaging these two with the two accurate ones will bias the corner position outward.
I'm thinking maybe the right procedure is not to average these values, but rather to trust the side that is closer to the center pixel. So the central pixel gets all of its own values from the original file that Craig verified. Then the 4 orthogonal pixels use the central pixel's values for the sides they share with it, but their own values for the rest. And so forth. Always trusting the values for the pixel with the smaller |i| or |j| as appropriate. I think this might give more accurate pixel areas. Especially (0,0) whose area should be literally the same.

[rmjarvis]

As a follow up to this, I think if we switch to a method where the 0,0 pixel is truly identical to the old code, we should add corresponding tests for (0,1), (1,0), and (1,1) (i.e. the next 3 most important pixels to get right). We should at least know how much their areas change between the old and new descriptions.

[jamesp-epcc]

This is now resolved after modifying the code to store extra corner points in the linear boundaries.

[jamesp-epcc]

Thanks Mike for the detailed review. On a first look through there's nothing there I really disagree with. I'll hopefully get onto addressing these points later in the week.

[jamesp-epcc]

I've pushed a new commit that addresses the simpler points from the review. The others will need a bit more work, especially the area issue. I will hopefully make some progress on them over the next few days.

[jamesp-epcc]

I've changed the initialization code as we discussed so that instead of averaging the points it uses the ones closest to the center pixel and ignores the rest (both for the linear and polygonal boundaries). This reduces the pixel area discrepancy from ~0.02 to ~0.005, but I don't really have a feel for whether this is close enough. I also had to modify some other expected results in the area test to make it pass, again I'm not sure whether these changes will be acceptable or not.

Strangely, I'm seeing test_flat fail on my system now (specifically the line that checks the variance), but it passes in all the GitHub CI checks, so I'm not sure what's going on with this.

[rmjarvis]

Why is the pixel area different at all now? The test was for the central pixel, which should have precisely the same boundary as it used to, no? I think some of the non-central pixels could have a slightly changed area, but I expected the central one to be spot on what we used to get. I think that's worth digging into further to try to understand.

[jamesp-epcc]

I am not sure it follows that the center pixel will necessarily have the same area just because the center distortions are the same as before. Because the pixels around the center also have a charge, there will be multiple distortion updates that affect the center pixel and most of these will be centered on other pixels, so the distortion values used on the center pixel will not always be from the center of the distortions.

[rmjarvis]

OK. That's fair. I guess a good sanity test though would be to add charge just to the center pixel. Then I think we expect the resulting area of the center pixel to be identical to what the old code gave, right?

And we can then look at the differences for the surrounding pixels to see how much they are. They should all be relatively small I think, but if any of them aren't we should make sure we understand why.

[jamesp-epcc]

It turns out there is already a test in test_silicon_area that has charge only in the center pixel, so I added some extra asserts to check the area of the center pixel and the surrounding ones. As expected, the center pixel area matches what the original code produces, but the ones surrounding it are slightly off. I've put the values produced by the original code in comments for comparison. The pixel areas to the left and right differ by about 0.001 while those above and below differ by about 0.0002. I don't really have a feel for whether that's close enough to be acceptable or not.

[jamesp-epcc]

I have removed the old code by wrapping in "#if 0" as discussed. Hopefully this is ready to be merged soon now.

[rmjarvis]

Code looks good James. Just two minor administrative items.

1. The file make_sensor_data_symmetrical.py is currently in the top-level directory. I think this belongs in devel or even devel/brighter-fatter maybe. This isn't a script we need to run on a regular basis or anything, right? It was just something you used to help confirm the accuracy of the code changes. So that's a devel kind of thing.
2. Can you please rebase this onto the current main branch? Something here became incompatible with the latest version of main, so it can't be rebased and merged anymore.

[jamesp-epcc]

I've moved the script into devel and merged the latest changes from main so that this should apply cleanly (I didn't do an actual rebase since I previously merged commits from main into this branch, but I can do if you prefer).

[rmjarvis]

Please rebase it. We try to avoid merge commits in GalSim unless absolutely necessary.

[jamesp-epcc]

No problem, have now rebased.

[rmjarvis]

Thanks James!