Weighted correlation matrix fixes

[The9Cat]

Summary

Computation of the weighted correlation matrix relies on the user to have computed a reconstruction with all PCs of interest. For example, if the last reconstruction used a particular pair of PCs, a subsequent call to wCorr() or wCorrAll() would be incomplete.

This addresses #114

Fix

• Calls to wCorr(), wCorrAll(), and wCorrKey() trigger a reconstruction for all PCs before computing the weighted correlation matrix
• These calls now take an optional nPC variable to limit the number of PCs used in the reconstruction. This defaults to the nPC value provided to the constructor
• To prevent recomputation on successive calls, the previous reconstruction will be used if the last reconstruction call was performed by a wCorr* call and the desired nPC is less than or equal to the number of the last call.

Checks so far

Checked using the pyEXP coefficients notebook.

[michael-petersen]

Following the example pyEXP_coefficients.ipynb, I do get the expected ssa.wCorrAll() behaviour, and if I put in an nPC value (i.e. ssa.wCorrAll(nPC)), I get the expected changes to the tutorial (note to self: change tutorial when this is complete!).

However, if I do

mat = ssa.wCorr("star disk",[2,0],12)-ssa.wCorr("star disk",[2,0],2)

I get all zeros -- I think something in the logic about recomputing is overriding the change. It may be that we want this behaviour (i.e. successive calls where one decreases nPC is a niche case), and one can override this by doing a fresh call to ssa.reconstruct(), but maybe there is some better solution?

Otherwise: looks good, works as expected, almost ready to merge!

[The9Cat]

Can you check the code box in the previous comment?

[The9Cat]

Okay, I can't tell what failure mode you are seeing that gives all zeros since your code block is blank. I tried it with various values of nPC and it seems to be as expected.

[The9Cat]

Oh, I think I've deduced what you meant: the code is still returning the degenerate part of the rank. Right, that's confusing. I'll check in a fix.

[michael-petersen]

Can you check the code box in the previous comment?

Bad formatting -- I've fixed it now. You fixed the issue in the code block, but no in the way that I was expecting! Either way, this seems sorted now, and I think the last to-do is figuring out how to make sure this is all clear in documentation.

[The9Cat]

Hm. In the example that you provided above:

mat = ssa.wCorr("star disk",[2,0],12)-ssa.wCorr("star disk",[2,0],2)

the current solution would fail because the ranks would be different. For your example to work, the rank of the matrix needs to be the same as the original npc at construction time with the matrix zeroed except for the active PCs specified by nPC. We could do that. However, that wouldn't be exactly what you want since the first wCorr call would include the weighted correlation of PCs 0, 1 with PCs >= 2. So your difference matrix would have wings. I think the user who typed your example wants to limit the number of PCs in the weighted correlation matrix from the top and the bottom.

We could revert to the previous form (with zeroing). Not sure...

[michael-petersen]

No, I don't think there is any need to revert. I think truncating the matrix is the right solution here, that way the user has a clearer example of what they are getting!

[The9Cat]

Okay. It's less flexible perhaps, but better defined. Specifically, let S(i) and S(j) be the series reconstructed from PC i, j and let (S(i), S(j))_w denote the Frobenius norm of the two trajectory matrices. Then, you get a nPC x nPC rank matrix with

wCorr(i, j) = (S(i), S(j))_w/( sqrt((S(i), S(i))_w)*sqrt((S(j), S(j))_w) )

[The9Cat]

I've documented the new behavior in the Python doc strings and defined the w-correlation matrix in the latest working version of the readthedocs manual.

Should we change the merge target to devel? Or leave it as main? I'm thinking it goes to main since the change fixes a behavior that might be called a bug.

[michael-petersen]

I agree that this can/should go in main; I've incremented the patch number to 7.8.1. I've been testing personally and the changes all behave as expected.

I've also used this as an opportunity to post the up-to-date PDF version of the paper (the markdown was already updated).