dstndstn / tractor

The Tractor: measuring astronomical sources via probabilistic inference

the Tractor

Probabilistic astronomical source detection & measurement

authors & license

Copyright 2011-2023 Dustin Lang (Perimeter Institute) & David W. Hogg (NYU/Flatiron)

Licensed under GPLv2; see LICENSE.

install

First, install the Astrometry.net code. (https://astrometry.net/downloads or https://github.com/dstndstn/astrometry.net/tags). You can do this by grabbing the code and using "make", or you can install the python code directly using pip:

pip install git+https://github.com/dstndstn/astrometry.net.git

Then grab the Tractor code:

git clone [email protected]:dstndstn/tractor.git
cd tractor

make
# OR
python3 setup.py build_ext --inplace --with-ceres --with-cython

It is possible to run directly out of the checked-out tractor directory. But if you want to install it, you can use pip, optionally with flags to enable Ceres Solver (requires the Ceres library), and Cython:

pip install -v --install-option="--with-ceres" --install-option="--with-cython" .

There is a test script that renders SDSS images:

python examples/tractor-sdss-synth.py

Prereqs:

• scipy
• numpy
• astrometry.net

Other packages used in various places include:

• fitsio: https://github.com/esheldon/fitsio
• emcee: http://dan.iel.fm/emcee/current

documentation:

Horribly incomplete, but online at http://thetractor.org/doc and http://thetractor.readthedocs.org/en/latest/

collaboration:

We are trying to move from a "Research code that only we use" phase to "Research code that other people can use". We are happy to hear your feedback; please feel free to file Issues. And naturally we will be pleased to see pull requests!

tractor's Issues

"system" paper abstract needs Lang/Hogg input

Well, the whole paper needs input. But first, the abstract!

Go through Tractor.optimize() carefully

eg, getChiImage() is called even though we may only need a subregion. If we have modelMasks, computing ROIs is easy.

directory structure

The top-level directory for tractor is a mess. I have no idea what to do there! But it seems to me that everything .py should probably be in a subdir and the thing should be packaged a bit.

installation failure because of missing mix_wrap.c

If I check out the package and try to install using the instructions given I get an error because of a missing file tractor/mix_wrap.c:

[djbard@ki-ls07 tractor]$ python setup.py install --prefix=$PWD/install/
running install
running build
running build_py
running build_ext
building 'tractor._mix' extension
C compiler: gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC

compile options: '-UNDEBUG -I/afs/slac/g/ki/software/python/2.5.4/lib/python2.5/site-packages/numpy/core/include -I/afs/slac/g/ki/software/python/2.5.4/include/python2.5 -c'
gcc: tractor/mix_wrap.c
gcc: tractor/mix_wrap.c: No such file or directory
gcc: no input files
gcc: tractor/mix_wrap.c: No such file or directory
gcc: no input files
error: Command "gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -UNDEBUG -I/afs/slac/g/ki/software/python/2.5.4/lib/python2.5/site-packages/numpy/core/include -I/afs/slac/g/ki/software/python/2.5.4/include/python2.5 -c tractor/mix_wrap.c -o build/temp.linux-x86_64-2.5/tractor/mix_wrap.o" failed with exit status 1

Indeed, tractor/mix_wrap.c is missing.

Try a DES Stripe82 brick

Point legacysurvey.org at the DECaLS site

Reduce "bands" to just the bands that actually have images.

Name our Great3 team

I proposed "Blamm!" (nominally Bard Lang and Meyers Marshall).

Support this or counter-propose!

Write coadd images out earlier

False Color Images

As I mentioned in the Inference group conference call, for debugging/paper purposes I've duplicated the SDSS's false-color image generation pipeline in Python. It can create images based on either a catalog or the loaded pixel data, if passed a Tractor as an argument. Here are examples from loaded pixel data and an inferred catalog from our MCMC work, respectively:

I need to clean up the code and comments, and fix it to work with Tractor instances that have more than 3 Images, but I'll post the code here in the meantime. If someone else wants to improve it or take a stab at it, feel free!

def falseColorImage(tractor, useModel=False):
  """
  Returns a PIL image which can be saved with rval.save(fname)
  or shown with matplotlib.pyplot.imshow(rval).

  Mimics the SDSS3 photoop pipeline stage sdss_frame_jpg.pro,
  which in turn calls djs_rgb_make.pro and polywarp_shift.pro.

  Outputs a false color image by mapping the i, r, and g bands
  to the RGB color channels, respectively.
  Because the images can have slight offsets due to the camera
  properties, the i and g bands are shifted so that they
  have the same centers as the r band.
  In the original pipeline, RGB channels are saturated at 30.
  (contiguous saturated pixels are replaced with their average 
  color), then scaled by 0.4 * [5., 6.5, 9.], then put through
  a nonlinear function [r,g,b] *= asinh(N*(r+g+b))/(N*(r+g+b))
  where N = 2.25.

  We simplify this by removing the nonlinearity and saturation.

  (Interesting tidbit from code: it seems there are 1361 unique
  rows in each image.)

  For now, we assume there are only images for the specified bands.
  """

  refimg = getFirstImgWithBand(tractor, 'r')

  # sdss_frame_jpg.pro
  # which does a polywarp shift of degree 1 (unless "keyword" set?)
  # filling with constant 0
  from scipy.ndimage.interpolation import shift

  def imgToCorrectedData(img):
    if useModel:
      data = tractor.getModelImage(img)
    else:
      data = img.getImage()

    return shiftBandImageToRef(img, refimg, data=data)

  idata = (0.4 * 5.) * imgToCorrectedData(getFirstImgWithBand(tractor, 'i'))
  rdata = (0.4 * 6.5) * imgToCorrectedData(refimg)
  gdata = (0.4 * 9.) * imgToCorrectedData(getFirstImgWithBand(tractor, 'g'))

  # print np.max(gdata)

  stacked = np.dstack((idata, rdata, gdata))

  stacked = np.flipud(stacked)

  # Perform saturation on individual pixels
  stacked[stacked < 0.] = 0.
  stacked[stacked > 30.] = 30.

  nonlin = 2.25
  sum_stacked = np.sum(stacked, axis=2)
  nonlin_mult = np.arcsinh(nonlin * sum_stacked) / (nonlin * sum_stacked + 1e-7)
  stacked *= nonlin_mult[:,:,np.newaxis]

  # Perform RGB saturation, which should be minimal
  stacked[stacked > 1.] = 1.

  from scipy.misc import toimage
  return toimage(stacked / np.max(stacked) * 256.0, channel_axis=2)

def getFirstImgWithBand(tractor, bandname):
  for img in tractor.getImages():
    if getBandNameForImage(img) == bandname:
      return img
  raise Exception("Tractor must have an image in band %s" % bandname)

def shiftBandImageToRef(img, refimg, data=None):
  if data is None:
    data = img.getImage()

  refpixctr = [n/2 for n in refimg.getShape()]
  refpos = refimg.getWcs().pixelToPosition(*refpixctr)

  from scipy.ndimage.interpolation import shift

  assert img.getShape() == refimg.getShape()

  pixdiff = img.getWcs().positionToPixel(refpos)
  pixshift = (refpixctr[1]-pixdiff[1], refpixctr[0]-pixdiff[0])

  shifted = shift(data, pixshift, order=1,
                  mode='constant', cval=0.0, prefilter=False)

  return shifted

Tractor Paper submission milestone needs a due date

Set up Travis-CI build

Installing Ceres might be a challenge...

Set up shared account with Tractor codebase

Get SLAC user account for Dstn

"Guards! Bring me the forms I need to fill in to have this person taken away"

• Futurama

@djbard, we need to send @dstndstn the link to the SLUO form and help him through that process. Would you mind please?

Read model-fitting shape estimation papers

Bridle et al started trying to measure shear by fitting blobs; Kuijken et al were not far behind. These days Bridle's method goes by the name Im3Shape and is calibrated to sims to solve the ellipticity prior problem. In between, Miller, Kitching et al applied some clever tricks to their own likelihood-based method, Lensfit, and applied it to the (multi-epoch) CFHTLens data. Lately, Bernstein has got interested in simple model fits, and their limitiations - with Armstrong. Links to follow, @drphilmarshall assigned. Bibtex file to be pushed before this can be closed.

problem running great3.py over a different image

I thought a good way for me to get more familiar with galsim and with the tractor would be to produce my own test images using galsim, and run the tractor over them. To this end, I made an image and a PSF (using demo2.py) and tried to run tractor on it. I'm running into problems with the model optimisation, and I could do with some help. The PSF and noise models look OK, but the optimisation fails entirely (it's 'borking'). Are there perhaps some defaults I should be tweaking? I've fiddled with everything I can see (initial params etc), and I've tried with different galaxy images (gaussian, sersic combos, high SNR, low SNR). I can't get anything to converge.

You can find an image and a psf file, the script to run tractor, and the logfile of the output from the tractor here:
http://www.slac.stanford.edu/~djbard/LSST/great3/
You should just need to edit the location of the files in great3-test.py and it will run for you.

README contains copyright statement

Looks inconsistent with the LICENSE file?
Could be nice if the README had some encouraging words for how to collaborate on tractor development and use?

Compare modelMasks to blob masks

ellipticity parameterization

We should think about (and implement!) different ways to parameterize ellipses, and think about which will be most robust.

Current parameterization of {radius, axis-ratio, position angle} is degenerate at axis-ratio = 1, for instance. Alternatives could include second moments {Ixx, Ixy, Iyy}, or radius/ellipticity {r, e1, e2} for instance.

Cannot run examples: ArithmeticParams not defined

Hi All,
I'm just trying a fresh install of tractor/astrometry, but when I try to run the examples I get this:

[djbard@ki-ls09 examples]$ python wldemo.py
Traceback (most recent call last):
File "wldemo.py", line 6, in
from tractor import *
File "/u/br/djbard/LSST/great3/tractor-master/tractor/install/lib/python2.7/site-packages/tractor/init.py", line 5, in
from motion import *
File "/u/br/djbard/LSST/great3/tractor-master/tractor/install/lib/python2.7/site-packages/tractor/motion.py", line 1, in
class Parallax(ArithmeticParams, ScalarParam):
NameError: name 'ArithmeticParams' is not defined

Am I missing something in my setup, or is this a bug?
-Debbie

Simplest Possible milestone has no date

What's a sensible date to aim for to get ourselves onto the leaderboard?
I guess it depends a lot on the SLAC set up time.

In tractor.Image, avoid storing both "invvar" and "inverr"

Choose one to be canonical and store only that. Probably "inverr".

Ensure that numerical derivatives have the same footprint

In some cases (maybe source radius set too small?; spatial derivatives) the footprints can be different and there can be significant flux at the edges, leading to large derivatives at the edges of the footprints

Set up multiprocessing CPU usage tracking

Infer the ellipticity distribution from the deep fields

Great3 includes deep versions of the fields. Use these to infer an ellipticity distribution (in preprocessing).

Test COSMOS region: bricks with tons of coverage

343924 149.12500 2.2500000
343928 150.12500 2.2500000
343929 150.37500 2.2500000
345366 149.62500 2.5000000
345371 150.87500 2.5000000
346806 149.62500 2.7500000
346807 149.87500 2.7500000
346810 150.62500 2.7500000

Final DR1 checklist

• Change "DECALSDR" FITS header to "DR1"
• Tag code
• Update output directory in pipebrick.sh

Hierarchical inference on the cheap

Suppose we assign default (ie uniform) Tractor priors to all parameters, including e1 and e2, and draw N posterior samples for each image. Can we implement a layered model after the fact to infer the intrinsic ellipticity distribution and shear field simultaneously? If we know the shear is constant in sky patches that will make life easier; if we can postpone the top layers of inference we can run in a pleasingly parallel way (and so do more experiments).

Make tests actually test things, or remove them

Match on David's photometric-cut CCD list.

/global/project/projectdirs/cosmo/work/decam/cats/ZeroPoints/ZeroPoints-DR1.fits

missing 'catalogs' in astrometry

When I compile astrometry, it fails due to a missing 'catalogs':

[djbard@ki-ls07 astrometry]$ make py
make -C qfits-an
make -C gsl-an
make[1]: Entering directory /nfs/slac/g/ki/ki06/lsst/djbard/LSST/great3/tractor/tractor/astrometry/gsl-an' make[1]: Nothing to be done forall'.
make[1]: Leaving directory /nfs/slac/g/ki/ki06/lsst/djbard/LSST/great3/tractor/tractor/astrometry/gsl-an' make[1]: Entering directory/nfs/slac/g/ki/ki06/lsst/djbard/LSST/great3/tractor/tractor/astrometry/qfits-an'
make[1]: Nothing to be done for all'. make[1]: Leaving directory/nfs/slac/g/ki/ki06/lsst/djbard/LSST/great3/tractor/tractor/astrometry/qfits-an'
make -C catalogs
make: *** catalogs: No such file or directory. Stop.
make: *** [py] Error 2
[djbard@ki-ls07 astrometry]$ ls
. .svn LICENSE PORTABILITY README.in init.pyc configure gsl-an make-snapshot.sh sdss
.. CREDITS Makefile README init.py astrometry freebsd-notes libkd qfits-an util

I got the astrometry package via the checkout.sh script.

Create general-Sersic index galaxy model

Probably fit the MoG approximations on a grid of Sersic n, and then interpolate MoG parameters.

SFD maps

We should probably divide them out at modeling time, if we can.

Outline the Blamm Great3 paper

We agreed to document what we're doing in paper form. Please create a .tex file and add sections for Great3 submissions we want to make. It can do in a new documents/ directory.

Run a test brick in our z-band-only region

Use ellipticity distribution as prior

Make the galaxy shape parameters use the (shared) ellipticity distribution as a prior.

Need to find a 'forum' to work in

We need some kind of working space that we can use to share updates on work-in-progress and discuss results. It needs to be able to host uploaded images etc, and allow comments and discussion. Organisation by tasks on a to-do list would be ideal.
Github isn't really the appropriate place place for this kind of detailed discussion. Options include: Trello, google+ circle, wiki page (hosted at slac?). Preferences? Other ideas?

Use CPDES, not DES

Contact Shreyas to set up a new Science Gateway for the interactive browser

getMag/setMag functions for NanoMaggies use different units

This is more of a "prepare this so it can be used by others" type issue, but still important enough to track IMO. In the process of tracking down some bugs for the MCMC project, I found that the design of the NanoMaggies class had a largely unintuitive and surprising gotcha; namely, a getter and setter use different units! Here's what I've found for how the tractor.basics.NanoMaggies class works:

• created in nanomaggies (linear)
• stored in nanomaggies (linear)
• print to strings in mags (log scale) UNLESS the nanomaggies happen to be negative, in which case it prints nanomaggies?
• exposed in nanomaggies (linear) by brightnessToCounts
• exposed in nanomaggies (linear) by getFlux and getBand
• exposed in mags (log scale) by getMags
• BUT updated in nanomaggies (linear) by setMags (defined in the superclass)

The relevant code is https://github.com/dstndstn/tractor/blob/master/tractor/basics.py#L206 - note that the setter from the superclass Mags is never overloaded.

For backwards compatibility with the rest of the code, we should be wary about changing this behavior, but since it is so unexpected, it should be documented much better. Alternately, if nobody is using setMags ever, we could add a setter that uses units of mags. I'll put together a pull request at my earliest opportunity, but I should collaborate with someone with more astronomical knowledge about the significance of these units before we change the documentation on master.

Feature request: FITS I/O of tractor Image objects

Will have to develop an API that allows each component of the image's calibration (PSF, WCS, Photocal, etc) to write its parameters to the FITS headers, and also a way of saying what type they were, and a way to recreate the object upon reading.

Add ln prior to tractor outputs

Now that tractor can apply an interim prior, can the ln prior be reported along with the samples?

Cannot install tractor!

I can no longer install tractor. Checking out a fresh version I can make it just fine but when I try to install I get:

[djbard@ki-ls09 tractor]$ python setup.py install --prefix=install
Traceback (most recent call last):
File "setup.py", line 12, in
eigen_inc = check_output('pkg-config --cflags eigen3')
NameError: name 'check_output' is not defined

Testing on GREAT08

Zuntz et al validated im3shape on the GREAT08 dataset, for which truth values are available. Seems sensible to experiment on this before submitting to GREAT3?

Using Astrometry without checkout.sh

I guess this is more an Astrometry bug, but since it's relevant to Tractor HEAD's use of Astrometry HEAD, I'm filing it here since Astrometry doesn't have an issue tracker. It shouldn't be necessary to create a local copy of Astrometry in a subdirectory the way checkout.sh does it. With a slight tweak, I was able to get Tractor (which I'm using as an installed package) to use the Astrometry package I installed through Homebrew according to the instructions here using --HEAD. Specifically, I needed to add util/ttime.py to the list of files that Astrometry.net installs in its makefile; it just needs to be added to the line in util/Makefile that includes yanny.py. In Homebrew syntax that's:

inreplace 'util/Makefile', 'yanny.py', 'yanny.py ttime.py'

So in other words, Astrometry added a new file that Tractor uses, but it was never added to Astrometry's build system configuration, so it doesn't get installed when Astrometry is installed.

Not a high-priority issue, but useful to track as we add collaborators.

Fix galaxy models going to crazy radius/ellipticity

Make sure to .gz appropriate image files

Fire up calibrations

Bring SLAC computing gurus up to speed

Tom Glanzmann and Joanne Bogard have agreed to help us run at scale at SLAC, with Richard Dubois overseeing. (This is a good example use case for us.) THey need to be brought up to date with the challenge and our plans/needs. @drphilmarshall to get this started, @djbard to meet them in person next week?

Add David's list of SDSS columns; and change output order.