Fix 105 and 118 #124
Merged
Fix 105 and 118 #124
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
As part of fixing #105.
"fitting.moments_enhanced" without "max_pix_variance_factor": removed from 1) guvectorize signature (two instances) 2) arguments 3) docstring 4) formula for "errorpeaksq"
6 instances of "max_pix_variance_factor" removed from image module.
"maximum_pixel_method_variance" no longer used, since its effect on the error analysis for peak spectral brightnesses is minimal.
"maximum_pixel_method_variance" and its return value, "max_pix_variance_factor", are no longer needed.
The "maximum_pixel_method_variance" method has been removed from PySE, so it does not need to be tested any longer.
1) Adds the "beam" (3-tuple) parameters to "fitting.moments" as an argument, to align it more with "fitting.moments_enhanced", which is considered superior for some edge cases. 2) Added "beam" to the docstring. 3) Fixes bug "maxpos = np.unravel_index(np.abs(data).argmax(), data.shape)". "np.abs" should not be there, could give wrong results for fits fixed to a position, i.e. without thresholding. 4) "pixels" docstring entry for "fitting.fit_gaussian" not completely correct. 5) "beam" docstring entry for "fitting.moments_enhanced" not completely correct.
1) Accommodate for extra argument "beam" in "fitting.moments" in the "test_gaussian" module. 2) "maj" and "min" are replaced by "smaj" and "smin" since we are dealing with the semi-major and semi-minor axes. 3) Adhere to a maximum of 80 characters per line.
1) The NVSS paper corrects for a general fitting bias, see equation (34) of the NVSS paper, as if fitted peaks are systematically biased too high. We don't observe that, rather the opposite, i.e. fitted peaks are biased too low, but more investigation is needed to quantify that. 2) The variable "rho3" is not used, so is removed. 3) Add comment that this assignment reflects equation (37) from the NVSS paper. 4) Since the peak spectral brightness is corrected, albeit from now on for clean bias only, this should be propagated in the flux densities. This is effectively equation (35) from the NVSS paper. 5) "data" in "extract.source_profile_and_errors" can be a "numpy.ma.MaskedArray". 6) "flux" should be equal to the addition of all pixel values from an island. 7) If some parameters are fixed, this should be reflected by the initial guess of these parameters. 8) A complete moments computation is only possible if we have imposed a non-zero threshold and no parameters are fixed.
When sources are fitted we compute the flux density from the peak spectral brightness and the fitted axes using the equivalent of equation (35) from the NVSS paper. Until now, we used the sum of all pixel values in the island when deriving an initial "guess" for the source parameters, but this too should be aligned with that equation (35). To be clear, moments are not yet computed in this part of the code, this is only an initial "guess" for the source parameters, which will be retained when moments computations fail. Effectively, this will mean that the flux density will be equal to the peak spectral brightness unless "beamsize" is computed in a "deviant" manner.
Removal of unnecessary conditionals in "fitting.moments" and "fitting.moments_enhanced". In order of appearance from a "git diff": 1) Clarification in docstring 2) "data.mean()" is probably not the right quantity in assessing whether we are fitting a positive or negative Gaussian, it can be flawed, e.g. when "data" covers a large patch of the sky, while "data.max()" cannot. 3) Wrt 2: "peak *= fudge_max_pix_factor" is appropriate in both cases. 4) There should be a check at some point in "fitting.moments" and "fitting.moments_enhanced" for the threshold to be below "data.min()" or we cannot correct the axes for bias from the threshold in an appropriate manner. The check from this commit does not work, i.e. we encounter "threshold >= data.min()" frequently. This is most likely because threshold is a single float equal to the analysis threshold map value at "maxpos". We will fix that in an upcoming commit. 5) The idea here is that the "if len(data.nonzero()[0]) > 1:" conditional is not needed since the case of a single pixel detection will be caught by either the "if semiminor_tmp <= 0:" conditional or the "if abs(semimajor_tmp - semiminor_tmp) < 0.01:" conditional. 6) If we have "semiminor_tmp <= 0" we can return the parameters equivalent to the "initial guess" immediately. The decision tree for calculating the position angle should be unaffected. 7) If we have "ratio < 1", i.e. "if threshold < peak", which should always be the case, "semimajor_tmp" and "semiminor_tmp" should always be positive and both "semimajor = math.sqrt(semimajor_tmp)" and "semiminor = math.sqrt(semiminor_tmp)' do not have to be encapsulated in a "try-except" block. 8) The computation of the flux density from moments should be aligned with its counterpart in Gaussian fits. 9) "target='parallel'" gives a speedup of just over 100 ms for 4K *4K images with ~167.000 sources. 10) Other changes in "fitting.moments_enhanced" are similar to those in "fitting.moments", with the exception of "min_with". In "extract.source_profile_and_errors" "fitting.moments" is called after checking for the minimal width. For its vectorized counterpart, i.e. for "fitting.moments_enhanced", that would not work, since the minimal width can be different for every detection, that is why this check is internalised. 11) As for "fitting.moments" the "threshold check" - now set at "assert threshold < maxi" or "assert threshold >= maxi"- should be adapted. It is now not a meaningful check while some threshold check is needed in order to correct for the axes, which are biased low from the threshold imposed for segmenting the island. "ratio = threshold/peak" in used to do so, it should be < 1 in any case.
It must have been at least a decade ago that a direct import of "ndimage" was possible, therefore the try..except block seems redundant.
Apparently, the update of this call to "fitting.moments' to accommodate for the extra "beam" argument had been missed. This is now fixed.
As we know, the axes from moments - when computed in the simplest way - are biased low. This is an effect from the threshold that segments the island. We can correct for that, but that requires all pixel values within the island to be above the (local) analysis threshold. That is what the "assert" statement in "fitting.moments" from a previous commit tried to accomplish. However, you could consider this a redundant check, except for edge cases where the island is selected not by segmentation, but e.g. by selecting a patch of the sky around some set of celestial coordinates. Some unit tests select chunks of observational image data in that manner and we have to accept that "fitting.moments" will return biased Gaussian axes in those cases. Still, the sanity check "peak > threshold" is required or we will end up with math domain errors.
1) The "extract.ParamSet().moments==True" condition to call "extract.ParamSet.bounds" is externalised, i.e. no longer internal to "extract.ParamSet.bounds", such that users that know what they're doing can call this method even if "extract.ParamSet().moments==False". 2) Always compute errors, either using "_condon_formulae" or "_error_bars_from_moments". 3) If any parameter is fixed, then do not compute moments. 4) Computing bounds - for Gaussian fits - only makes sense if no parameters are fixed, since imposing any bounds on non-fixed parameters will, most likely, result in a fitted value equal to one of the bounds, i.e. a compromised fit. 5) The case of "if fixed and not param.gaussian:" can be left out, since we can still proceed with the values from the "initial guess".
1) moments needed to compute bounds. Until recent, that was already the case, but that condition is now externalised. Also it was felt that, for the regular source measurement process, i.e. when "extract.source_profile_and_errors" is called, having no fixed parameters was a too weak condition; moments have to be computed succesfully in order to derive bounds in a meaningful way. 2) Clarification of comment.
"known_result_moment" had not yet been adapted for the removal of "max_pix_variance_factor". This has a tiny effect on the error bars on both the peak spectral brightness and the flux density.
The removal of the bias correction from "extract.Paramset._condon_formulae" affects the ground truth peak spectral brightness, its error bar, chi-squared and reduced chi_squared for the "testSingleSourceExtraction" unit test. This commit adapts that ground truth data to our code. Some background: we don't observe the "general fitting bias pulling the fitted peak in the direction of the local noise gradient" as reported in the NVSS paper, equation (34), at least not as a positive bias. We have therefore removed the second term on the right-hand side in the line of code that mimics equation (34), in a previous commit. The new ground truth data shows both a lower chi-squared as well as a lower reduced chi-squared.
The ground truth data for the peak spectral brightness had been set at -0.00364 Jy/beam. It turned out this number came from negative bounds from a moments analysis; the "fitting.moments" code used to derive that number still used "data.mean()" in the conditional for deciding to derive a negative or a positive peak. Since "data.mean()" was negative, the bounds were negative. In our updated code, we use "data[maxpos]" instead. Moreover, we do not compute moments when one or more Gaussian parameters are fixed, consequently no bounds are set for the "testSourceAtGivenPosition_negative_spectral_brightness" unit test; the peak fit is unbounded, which yields this new ground truth peak spectral brightness of -0.00053898 Jy/beam. This includes the removal of the bias correction, i.e. the removal of the second term on the right-hand side of our implementation of equation (34) from the NVSS paper.
Better use data.max() instead of data.mean() in deciding to derive either a positive or a negative peak. This mimicks the conditionals in "fitting.moments" and "fitting.moments_enhanced".
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Fixes #105
Fixes #118