Angs

[garrettwrong]

Stashing debug work.

So far I have found two errors in the initial estimating rotations method preventing us from achieving the results from MATLAB. First we should be taking the average of 1d (alpha) angles, not the rotation matrices. Second, the euler to rotation conversion in our current python code is not equivalent to the one in MATLAB (and originally ported by Junchao).

I saved the alpha angles from MATLAB for all i,j image pairs. Following #1171 we can make an equivalent CL mat. Following above we can take the alpha MATLAB angles and generate the same rotations.

There is a small (1.5*) discrepancy between the alpha angles computed in Python and those in MATLAB. This is coincidentally the value I would expect in this particular problem if we had an 'off by one bin' error. (Not sure, just thinking atm).

I will continue working this, hopefully resulting in MATLAB parity up to initial rotation estimates for one of the JSB problems.

[garrettwrong]

These changes are the minimal set that allow me to reproduce the initial rotations estimate from the published MATLAB common line sync 3n algorithm (up to transposing the image data). For the amount of effort to get there, its not much code to look at :D.

I have put in branches/options with different defaults to allow the other python common lines implementations to continue functioning as they were.

I was able to salvage using the ZYZ convention so we don't need to worry about that ATM. yay 🥳 .

[j-c-c]

Looks great! Can't believe it only took ~100 lines.

Just one small docstring thing. Also, looks like CI failed. Is that from flakiness or is something wrong there?

[garrettwrong]

To be clear, this is just to fix the voting procedure angle outputs. Gonna go on a limb and guess the rest of the CL algorithm(s) in the current code base are also not reproducing MATLAB.... one piece at a time...

Saw that unit test fail for first time last night. Not sure about that yet.

Because it is almost exactly 90* off, im wondering if it relates to the synchronization and convention I tried to put back. I got out of using the old rotation convention/conversion because the sync condition was "fixing" the cases that were different into the domain that is equivalent (but maybe the sync check it isn't working quite right numerically on that one platform). I'll try to look into that more.

[garrettwrong]

Yeah, I'm not really sure what is going on with the test. I can reproduce the failure locally.

Its almost like the remaining common lines code doesn't work with the correct initial rotation estimates.

[garrettwrong]

Yeah, I'm not really sure what is going on with the test. I can reproduce the failure locally.

Its almost like the remaining common lines code doesn't work with the correct initial rotation estimates.

Oof, apparently the upstream code doesn't like the patch that returns actual rotation matrices. Reverting allows the test to pass. I figured there was more to do here, but was hoping I wouldn't need to do in this PR.

diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index 8802ec05..7ee5cbad 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -815,14 +815,14 @@ class CLSync3N(CLOrient3D, SyncVotingMixin):
         """
         alphas, good_k = self._vote_ij(clmatrix, n_theta, i, j, k_list, sync=True)
 
-        angles = np.zeros(3)
+        angles = np.zeros((len(alphas), 3))
 
         if alphas is not None:
-            angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
-            angles[1] = np.mean(alphas)
-            angles[2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
+            angles[:,0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2
+            angles[:,1] = alphas
+            angles[:,2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta
             rot = Rotation.from_euler(angles).matrices
-
+            rot = np.mean(rot, axis=0)

Note the direct averaging of rotation matrices is used by every other CL algo we have in Python...

[garrettwrong]

Okay, had to dig in more and found an additional issue with nearest_rotations compared to MATLAB which I've patched here for Sync3N. I think patching that has fixed the unit test issue.

During this effort, I confirmed for a small (unit test size of 41px) problem that this code can now reproduce MATLAB sync 3N's output rotations up to global reflections when given an identical common lines matrix. This was done for cases S0J0 S1J0 S0J1 S1J1 corresponding to S_weighting and J_weighting. I will be moving onto the next problem size, 80s at 89px. Afterwards I will put this up for review consideration again.

[garrettwrong]

Reminder to squash this PR during merge, lots of junk debugging commits. squashed

new changes (nearest_rotation) seems to resolve unit test

[garrettwrong]

@j-c-c , can you take another look at this before I send it to Joakim, thanks!

[j-c-c]

@j-c-c , can you take another look at this before I send it to Joakim, thanks!

Sure!

[j-c-c]

Looks great!

[janden]

Looks great! Just one thing.

[janden]

Other way around, no? (Setting the option does not correct reflections.)

[garrettwrong]

Setting the option to True allows reflections.

I'll add disable to the docstring to better indicate the default; which I think is what you meant by setting.

[janden]

Thank you!