astrorama / SourceXtractorPlusPlus

SourceXtractor++, the next generation SExtractor
https://astrorama.github.io/SourceXtractorPlusPlus/
GNU Lesser General Public License v3.0
72 stars 9 forks source link

Save catalog even if model fitting does not handle all detections #334

Open andalenavals opened 3 years ago

andalenavals commented 3 years ago

Hello,

I am running sourcextractor++ on some simulations using FlexibleModelFitting and a constant psf (I have tried both fits file and psf-fwhm and psf-pixel-sampling values ). I am getting for several simulations that not all the detected objects are measured.

I have save all debugging levels (log-level) but It seems there is not a major explanation of the kind of error during the measurement. And there is not a catalog production for the objects that were successfully measured.

I am using https://anaconda.org/astrorama/sourcextractor 0.14 version in a conda env with python 3.6

Questions are, what is it the best way for tracking potential bad measurements beyond (log-level). And the not production of the catalog even if all the objects were not detected seems like a bug.

I am attaching the debug output debug.txt

ayllon commented 3 years ago

Hello @andalenavals

This looks like some kind of crash. You should get a line like X sources detected, or NO SOURCES DETECTED at the very least. However, the prefetcher and the output threads seem to finish orderly.

I would suggest to run the detection without model fitting, purely detection (asking only for PixelCentroid, for instance) This should help reducing the surface area, and let us see whether we have a detection problem, or a fitting problem.

Of course, whatever the error it should be either fixed or dealt with more gracefully.

ayllon commented 3 years ago

Hello again,

I think I have hit this issue while chasing something else. It seems that there is a race condition between when some dependent parameters raise a python exception, and when sourcextractor++ translates it into a C++ one. If you get unlucky, the raised error gets unset in the middle, and it gets lost. However, the execution is aborted all the same.

As a workaround, running with thread-count=1 may give you the error that is causing the stopping of the process.

andalenavals commented 3 years ago

Hello @ayllon , sorry I did not attach the INFO log, info.txt

here it is where the remaining output is saved, including how many sources were detected and measured (FlexiModel apparently measured 2 objects, but there is not output catalog of those). I have already set thread-count=1, in the debug file it appears as defined. But actually I see in my machine that more than one core is activated, I tried the trick of issue #333 but multiple cores are being called yet.

When I switch to PixelCentroid everything works fine, and I get the catalog with the pixel centroid of the objects, so it is not a detection problem.

Something else that happens after running with FM is that the terminal fomating is lost, I am unable to see what I am typing and new lines are not aligned

ayllon commented 3 years ago

It still seems like a crash to me. It shows 2/9, which means actually 9 where detected, and it probably crashed while fitting the third. This also explains why the terminal is unusable (btw, you can try typing reset).

Would it be possible to get a copy of your config files and the image you are fitting? The fact that at least 2 are measured discards the possibility of an exception due to a typo or the like.

andalenavals commented 3 years ago

Hello @ayllon , I have done a simplified version using constant psf.

standalone_Error.tar.gz

ayllon commented 3 years ago

Great, thanks! I can reproduce immediately the crash with the develop branch. I reopen, btw, I guess you didn't intend to close?

ayllon commented 3 years ago

And as I suspected, #336 fixes it. If I run with that branch, I get the error:

2021-04-15T13:08:50CEST Config ERROR : Python exception Latitude angle(s) must be within -90 deg <= angle <= 90 deg, got 98.03376770019531 deg
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 1034, in wc_rad_func
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 904, in get_separation_angle
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 868, in get_sky_coord
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate.py", line 255, in __init__
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 245, in _parse_coordinate_data
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 591, in _get_representation_attrs
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/angles.py", line 507, in __new__
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/angles.py", line 528, in _validate_angles
2021-04-15T13:08:50CEST Config ERROR : Python exception Latitude angle(s) must be within -90 deg <= angle <= 90 deg, got 98.11773681640625 deg
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 1034, in wc_rad_func
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 904, in get_separation_angle
2021-04-15T13:08:50CEST Config ERROR : File "/home/aalvarez/Work/Projects/SourceXtractorPlusPlus/SourceXtractorPlusPlus/SEImplementation/python/sourcextractor/config/model_fitting.py", line 868, in get_sky_coord
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate.py", line 255, in __init__
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 245, in _parse_coordinate_data
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 591, in _get_representation_attrs
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/angles.py", line 507, in __new__
2021-04-15T13:08:50CEST Config ERROR : File "/usr/lib64/python3.9/site-packages/astropy/coordinates/angles.py", line 528, in _validate_angles
2021-04-15T13:08:51CEST Multithreading DEBUG : Stopping output thread
2021-04-15T13:08:51CEST SourceXtractor FATAL : Latitude angle(s) must be within -90 deg <= angle <= 90 deg, got 98.11773681640625 deg

So at least that solves the crash without further explanation.

As for the error itself, it seems we assume sky coordinates to always be in degrees, and for this particular image it is not the case, it is just a direct translation pixel to pixel. We should probably check the CTYPE headers.

As a workaround, you can comment/remove these lines:

    ra, dec, wc_rad, wc_angle, wc_ratio = get_world_parameters(x, y, rad, angle, ratio)
...
    add_output_column('mf_ra', ra)
    add_output_column('mf_dec', dec)

Doing so it runs fine.

model

andalenavals commented 3 years ago

Thank you @ayllon. There is not crash now, and the catalog is being produced. However, running time is yet long. It took around 20 minutes in a machine with 64 cores to finish one simulation with 4096 objects. Any advice for this last point?

ayllon commented 3 years ago

Any advice for this last point?

@mkuemmel has found very recently that this might be due to a "conflict" between conda's libblas and sourcextractor++ multi-threading. They both try to use all cores, which adds too much overhead.

Try setting these variables on your environment and let us know:

export MKL_NUM_THREADS=1
export MKL_DYNAMIC="FALSE"
export OMP_NUM_THREADS=1
export OMP_DYNAMIC="FALSE"
mkuemmel commented 3 years ago

@andalenavals Yep, looks like you run on conda, and there it is really worth setting those variables. That could reduce the run times by quite a factor! Also using v0.14 would help (but I think you are already there).

We are currently discussing to set those variables within SE++ for the next version then.

andalenavals commented 3 years ago

Excelent!! What a huge difference. Setting those variables time reduced from 20 min to only 2 (using a single core). And when using 64 it took only 33 sec! Great. Thanks

mkuemmel commented 3 years ago

Looks like you have to still throw in some objects, otherwise your 64 cores get bored.

andalenavals commented 3 years ago

Hello again,

Let me continue a related issue I just found.

If you take the exact image I shared and configuration (after commenting the problematic lines). Then the produced catalogs give mf_rad of 27, 11 for a couple of galaxies, which by simply eye inspection you can see units are not in pixels. So I am wondering what units are being used here, is there any way to get values in pixels? Actually, there is even a more problematic issue and it is that comparing with true values I am getting a flat relation for some quantities, meaning the sersic index and the radius. Let me share you some plots explaining this last statement quicktestplot11 quicktestplot13

marcschefer commented 3 years ago

The default unit for the radius is indeed pixels but those are pixels on the detection image. We can handle multiple frames so the detection image grid is used as a common reference.

andalenavals commented 3 years ago

If that is the case, then I am getting not meaningful radius as you can see in plots above and also in my detection image, actually all mf_sersic_rad are above 16. Let me share my output catalog for those 9 sources in here

Hello @ayllon , I have done a simplified version using constant psf.

standalone_Error.tar.gz

cat.zip

ayllon commented 3 years ago

I think there is one configuration issue: if you combine Exponential and Sersic models, you need to give each its own flux parameter, since they will be both contributing to the total.

i.e.

    flux_total = get_flux_parameter()
    sersic_frac = FreeParameter(0.5, Range((0.0,1.0), RangeType.LINEAR))
    flux_sersic = DependentParameter(lambda f, r: f*r, flux_total, sersic_frac)
    flux_exp = DependentParameter(lambda f, r: f*(1-r), flux_total, sersic_frac)

    # We want the magnitude too
    mag = DependentParameter(lambda f: -2.5 * np.log10(f) + args.mag_zeropoint, flux_total)

    # Models to be fitted
    add_model(group, SersicModel(x, y, flux_sersic, rad_sersic, ratio_sersic, angle_sersic, sersic))
    add_model(group, ExponentialModel(x, y, flux_exp, rad, ratio, angle))

This does not seem to solve the radius, though. I need to double check with the others, but I think this might be related to #215 / #230. Up to a point part of the raster will be outside of the cutout used for the fit, so the radius may grow past it.

ayllon commented 3 years ago

I should have added: it may grow past it, and what lies outside the cutout will not contribute to the residuals.

ayllon commented 3 years ago

I tried adding a very simple prior and it constrains the radius:

    add_prior(rad, lambda o: o.get_radius(), 1)
    add_prior(rad_sersic, lambda o: o.get_radius(), 1)

rad_prior

The colored line corresponds to the Sersic profile (with the color being the Sersic index), the blue to the Exponential.

I guess the solution can go somewhere along that line.

andalenavals commented 3 years ago

@ayllon thanks. Now it looks better for that image, But let me add here more challenging ones.

img.zip

Could you try to repeat those plots, specially this one 20210519T214650_5ppgcsq7_img4_galimg.fits I believe in those you will observe more clear the effect I am talking about. Actually the plots of the scatter of sersic index and radius as a function of true values did not change much (with the updates you suggessted).

ayllon commented 3 years ago

It does work for me. No prior:

no_prior

With prior:

with_prior

This is the command I used:

sourcextractor++ --config-file constant_psf.conf --progress-bar-disable --output-catalog-filename=img4_cat.fits --detection-image 20210519T214650_5ppgcsq7_img4_galimg.fits image_file=20210519T214650_5ppgcsq7_img4_galimg.fits

Well, with "works" I mean I get no nonsensical radius. Might still have a weak correlation with truth, but that I can not tell.

andalenavals commented 3 years ago

I was able to reproduce the plots, I was pointing to the wrong python-config-file without priors. And as you suspect I am getting weak correlation plots yet. quicktestplot11 quicktestplot13 quicktestplot14

ayllon commented 3 years ago

Hello,

I can see another couple of things that may help:

  1. I do not think there is a background? It should be disabled if that is the case: background-value=0.0 on the *.conf file
  2. The ratio for both Sersic and Exponential should be in the range (0,1] (note that 0 is not in the range). Since it is an ellipse, you want the ratio between the major and the minor axes. If you let the value go beyond 1, you are making the fitting more difficult, since, at the very least, there will be two equivalent solutions (i.e. ratio 2 and angle 90 vs ratio 0.5 and angle 0)
  3. I think the ratio works better when using exponential. Honestly I do not know why, but that's what sextractor does too
  4. Maybe restrict the angle to a semicircle between pi/2 and -pi/2?

Example for the ratio:

    ratio = FreeParameter(0.7, Range((0.01, 1), RangeType.EXPONENTIAL))

You may have to play a bit with the initial conditions, ranges and priors. Ultimately, the fitting has no concept of right or wrong, only lower or higher residual.

residuals

residuals_hist

Those are the residuals with the parameters I have suggested. I can hardly blame the engine for thinking that's basically noise.

andalenavals commented 3 years ago

Great @ayllon, yes those images have no background I updated my configuration file. I also play a little with the limits and priors for the python-config-file (attached). Might be this tunning can be tricky but so far not an improvement on my attempt of getting correlated sersic indexes and radius with the true values. Might be those are noisy and hard to determine given the conditions. Anyways at least for the flux and snrratio things look better.

fitting_constpsf.py.txt quicktestplot00 quicktestplot05 quicktestplot13 quicktestplot11