darkenergysurvey / ugali Goto Github PK

This is an enhancement split out of #15. One question is how we want to deal with unbounded parameters. Adding @apace7 due to past interest in this topic.

@apace7 noted that the points we are using to seed the mcmc are not explicitly constrained to the allowed parameter range. This is a bug and should be fixed. @apace7, any chance you could submit a pull request with the fix you've been using?

FNAL, Midway, SLAC config files (better to make batch-system independent)

We could add a more complete and modern galaxy sample from HyperLeda. This is a continuation of the RC3 catalog of bright galaxies. I haven't noticed any major issues with incompleteness in the catalogs we have ingested (NGC, IC, and UGC), but I am continually paranoid about missing something.

Currently we only have Chabrier, we should at least add Kroupa.

Check for numerical issues with very small Plummer radii, sampling kernel (at the very least add a warning). Investigate large kernel fitting accuracy, what is largest size?

@bechtol has been developing simulation tools under the scratch directory. Once development slows it would be good to re-integrate these into the simulation module.

Now that we are using ugali across surveys it probably makes sense to start moving the package to a github organization. This would also allow us to host dependent tools such as the simple_binner, dsphsim, and a future "pysochrone" module.

It would be nice to add the background density to the results file. This is a bit complicated, because we are actually calculating the background model as a 2D probability distribution in color magnitude space, but maybe we could do something like sum loglike.b values for stars in the interior region?

This feature was requested by @scantu and @sidneymau in order to better plot radial profiles including the background density.

Residual maps / how to deal with gradients in background (important mainly for LMC region)

Distributing the isochrones is a headache. Might we want to compress isochrone data to easy distribution and allow for more extensive isochrone set? This is connected to #21.

Working on MCMC analysis with ugali and tyring to use g-i bands has shown that ugali.analysis.source seems to create a g-r isochrone even when the relevant config files specify band_1 = g and band_2 = i. This behavior has been evidenced in ugali/pipeline/run_05.0_followup.py and ugali/analysis/mcmc.py, but may be elsewhere.

Make mag_1_field more consistent (propogate 'g','r','i' throughout)

Be careful with spatialBin function when objects fall outside of ROI. This may actually be throwing an error now (which is not terrible). Need to look into more deeply.

We changed the default coordinate system to equatorial, but this didn't fully propagate to the results file. In particular, we are mis-reporting the glon,glat values and the iau name. Taking into account the correlated uncertainties between ra,dec and glon,glat is most easily done at the level of estimating parameters from the chains.

I don't know what this means anymore.

This is a moving target since new isochrones come out all the time. It's also a pain because of the way we distribute the isochrone files. We might want to think about a standalone isochrone package (certainly a lot of work to be done there). We might also want to think about a deeper integration of the code to download isochrones from the web.

@willcerny tracking this as a separate issue since #78 is closed. It should be pretty easy to calculated surface brightness from the Martin mag if it exists and output it as an additional results field. Basically add something like this inside the if martin block:

Precompute HEALPix pixel for each object and avoid creating full HEALPix maps. I'm not really sure what this means any more, but I think it would involve some pre-processing before catalog creation.

The column for log_lum should be 6 (as opposed to 8) for DECam based on an example isochrone from MIST.

The line to change is here:

Columns from the isochrone file:

#   1                               2                               3                               4                               5                               6                               7                               8                   9                  10                  11                  12                  13                  14                  15                  16
# EEP          log10_isochrone_age_yr                    initial_mass                       star_mass                        log_Teff                           log_g                           log_L                     [Fe/H]_init              [Fe/H]             DECam_u             DECam_g             DECam_r             DECam_i             DECam_z             DECam_Y               phase

Following on #34, the simple binner now lives happily in sidneymau/simple and can be removed from scratch.

In Line 1219 and 1220 of ugali.utils.plotting.py, an error is returned in calls to plotMembership because of a reference to self.isochrone.band_1 and self.isochrone.band_2. The "self" should be removed in both cases.

Also, the axes labels are swapped - line 1219 should be set_xlabel and 1220 should be set_ylabel

Matplotlib b2.2 has started warning that mlab.rec_append_fields is deprecated. There is a drop-in replacement in numpy.lib.recfunctions; however, it is a factor of 5x slower...

We need a more robust way to check for strings. I think that putting an isstr function in utils.fileio.py like suggested in this answer would be a reasonable solution.

We made the unfortunate (though physically motivated) decision to use Galactic coordinates as the default. We need to change this to Celestial, and make sure that the infrastructure for changing coordinate systems actually works.

This is a great space saver, but uses a non-standard file format. We should investigate whether healpy has evolved to have a sparse file format.

We are currently inconsistently defining whether our histograms are defined with hist.shape == (len(xbins),len(ybins)) or hist.shape == (len(xbins)-1,len(ybins)-1). The second defines the bin value at the center of the bin and is the convention followed by numpy.

It would be useful to define a catalog for diffuse structures (streams & clouds) that can be potentially confused with satellite candidates. I think the preliminary list would consist of:

• ATLAS stream (particularly the over dense "progenitor")
• Phoenix stream (particularly around over densities)
• SMCNOD
• LMC tidal tail from Mackey

We are sometimes seeing that the peak interval does not contain the peak (kde peak)
I think this has to do with lines 99 and 100, which currently read:
lo = x[min_idx]
hi = x[min_idx+window]

and should probably be changed to something along the lines of
lo = starts[select][min_idx]
hi = ends[select][min_idx]

We currently track two coordinate systems: Galactic (gal) and celestial (cel). Celestial just means "sky" and doesn't actually contain information about coordinate frame. The actual system is the equatorial coordinate system, so we should update instances of cel to equ. This is a trivial change, but permeates the codebase.

Investigate integration with astropy. Of specific interest is the astropy coordinate module; however, coordinate transforms are notoriously slow (astropy/astropy#1717)

The isochrone download scripts have become out of date as we've moved the server query functionality into the isochrone classes. We should update these scripts to use the current configuration in the isochrone classes.

Pinging @mcnanna to keep an eye on this for his LSST work.

Merge config file with default config or store defaults in objects (probably the latter)?

@sidneymau @bechtol do we want to pull out the simple binner into it's own repo that depends on ugali? Maybe @sidneymau did this already, or maybe it's more work than we want?

@sidneymau ran into the problem that the absolute magnitude calculation requires that the isochrone uses g and r band magnitudes assuming the DES system. We can work around either by having transformation equations for other bands, or by accessing the g and r magnitudes directly from the isochrone. The use of sample makes this a bit more complicated, so maybe we need to make a copy of the isochrone with band_1 = g and band_2 = r.

I've found fitsio to be superior to pyfits. In addition, astropy has taken over pyfits, which is now accessed as astropy.io.fits. I'd like to switch to fitsio, or abstract the interface to the fits reader so the user can use whatever is available on their system (fitsio > astropy.io.fits > pyfits).

Arguably, this line should be using the variable rsize (the azimuthally averaged half-light radius in arcmin) for the calculation (instead of size), since the SB formula is assuming a circular area.

The Marigo2017 isochrone tarball has two points which are incorrect. One looks to be the last data point in the file which has a discontinuity to the other points.

test.pdf

We implemented the ability to output the Galactic coordinates with the "correct" coordinate uncertainties by directly converting the mcmc chains. However, this never propagated through to the results file. In fact, the gal values in the results file have been incorrect.

Generalize bands in plot output, i.e.

Examples:

Pythonic style string formatting

plt.xlabel(r'${} - {}$'.format(self.isochrone.band_1, self.isochrone.band_2)
plt.ylabel(r'${}$'.format(self.isochrone.band_1))

C style string formatting

plt.xlabel(r'$%s - %s$' % (self.isochrone.band_1, self.isochrone.band_2))
plt.ylabel(r'$%s$' % self.isochrone.band_1)

Hahaha, what a joke... but the title says it all. Of course, this issue will never be fully resolved, but we can make a start.

The isochrone constructor currently disregards user-specified bounds on age/metallicity and sets them to the range of the isochrone grid (see here). Possible improvements:

1. Continue this behavior but warn the user
2. Respect the user's bounds when they are within the isochrone grid; when outside warn the user and set to the isochrone grid bounds.
3. Respect the user's bounds and throw and exception when they are outside the allowed grid bounds.

For both 2. and 3. above, the default isochrone grid consists of age = [10, 12] and z = [0.0001,0.0002], so the default bounds should be similarly constrained.

There is a fragility in the connection between the selection applied to the isochrones and the selection applied to the data. We want to connect isochrone.observableFraction and mask.restrictCatalogToObservableSpace better.

Compute observable fraction including photometric errors and detection completeness. Currently only using the magnitude limit. Fully incorporate completeness in isochrone fitting and simulation.

A few issues with the catalogs:
Laevens 1 is listed in the ExtraClusters catalog, but the ra, dec, and reference are all for Laevnes 2/Triangulum II. Triangulum II already exists in the ExtraDwarfs catalog. I think the solution is to change the ra, dec, and reference for Laevens 1 to actually be for Laevens 1.

Aquarius II is not in any catalog. It probably belongs in ExtraDwarfs.

McConnachie15 has an object named "Indus I ?", which is actually Indus II.

This may already be done in isochrone.py.

Ugh, but LSST is moving this way so we better too.

Calculating the per pixel surface intensity for small kernels has been a continual issue. The basic problem is that when the pixel size becomes comparable to the kernel size, the location of the kernel centroid within a pixel can have a large effect on the surface intensity calculated within that pixel. This leads to an unnormalized surface intensity and a mis-estimation of the richness and other source parameters.

One possible solution would be to subsample the surface intensity within each pixel (or just the central pixels). Whether this is technically and computationally feasible needs to be investigated.

This should be an easy one...