Would be possible to speed up the calculation more in the case of the SFH, IMF, and age-metallicity is fixed. Do not need to calculate the full physics model grid - only the points that are the in the sparse likelilhoods for the stars being fit. The code used to work this way, but this was removed when adding the fitting for the stellar parameters.

Set to by python 3.5 and above.
Remove earlier versions of python from travis.

Really need to add some basic tests a minimum. Currently there are no tests.

The core megabeast computation for a single ensemble run (e.g., single pixel) would be better in a callable function instead of where it is right now.

Current model limited to A(V) only.

Should consider finding off fields for each galaxy. If they don't exist in the archive (via some deep field), we should consider an observing proposal. Such off fields should be added to future proposals.

The BEAST has a fairly well-defined directory structure - from the folder it's run in, the data files are assumed to be in ./data/ (and new data-related files go there too), and the various modeling outputs are put in ./projectname/. The last BEAST step of doing the spatial reorganization lets the user put the files anywhere, and megabeast gets run from inside that folder. I'm thinking it would be nice to standardize that part too.

Uses BEAST physics and observation model.
Priors on ensemble (meta) model -> based on BEAST prior model.
Sparse likelihoods from BEAST runs.
Everything on this grid -> speed.

• Distance priors need to be flexible enough in MegaBEAST to allow for distance changes of the order a few tenths of a mag (maybe 0.5? more?) in distance modulus.

• Depending on how distances are fit in BEAST, the prior may simply be an updated best distance value. But, if BEAST is fitting distances based on specific features in a CMD, the input may be an updated zeropoint calibration. For example, if the horizontal branch was used in BEAST, then the zeropoint could be changed (or the uncertainties on the zeropoint). If the TRGB was used, the zeropoint could similarly be changed. If both features are used, then this is more difficult, and may require a re-run of the BEAST?

Inspired by @kmcquinn's ISM*-at-ST talk, the different parts of the CMD (HR diagram) are known to be better/worse understood (e.g., Hidalgo et al. 2011). This could be easily reflected in the megabeast fits by adding a model uncertainty term to the fitting. Or maybe this is a straight weight in the Bayesian formulation. Maybe @eteq can provide some advice here.

The completeness is a function of source density.
The current implementation only allows for a single noise model, hence source density. This should be generalized. Either by using the stats file to get the reorder_tag (I think) to say which noise model to use, or get add this info to the lnp file when reordering the BEAST results.

(this is an issue to put down ideas - useful before and during the actual coding)

To test the megabeast, need to have simulations that fully test the different parameters of the ensemble model. The simulations are for megabeast ensemble parameters (e.g., IMF slope, Av distribution). There is code currently to simulate A(V) ensemble models and it needs to be updated/generalized.

To start with, the simulations should be done with easy to represent uncertainties on the beast parameters with the beast pPDFs simulated from the physicsgrid. Later, the simulations can be made more realistic by simulating stars with fluxes and uncertainties using the beast code to do this simulation. The latter simulation is already implemented in the beast. Maybe it would be better to just do the more realistic simulations now. That's where we need to be in the end anyway.

One idea it to use the prior model in the beast as this is the ensemble model for the megabeast. The megabeast ensemble model can be specified using the same dictionary structure as the beast prior model.

The megabeast simulation code would need to (for a single pixel):

1. calculate the beast prior model based on the megabeast ensemble model based on a beast physics modelgrid
2. create a new physicsgrid with the new weights based on the ensemble model
3. use the beast simulation observations code to generate a simulated observed catalog
4. run the beast on the simulated catalog and the original beast physicsgrid model for a specified number of stars
5. run the megabeast on the beast results and solve for the ensemble model
6. compare the megabeast solution to the input megabeast ensemble model

Simulations are currently done in the BEAST code. This requires generating a new physicsmodel for each change in the BEAST priors.

The MegaBEAST work has provided code to update the physicsmodel "prior_weight" without needed to regenerate the full physicsmodel file. Moving the simulations to the MegaBEAST would speed up simulations. This update should use the megabeast settings format to specify the ensemble parameters (=BEAST priors) and update the prior_weight on a previously computed BEAST physicsmodel.

The full megabeast ensemble model needs to be implemented and documented. This includes allowing for the ensemble parameters to be input, the "simple" model implemented, and the details put in the docs.
Components of megabeast model:

• A(V) distribution (log-normal - confirm fully implemented)
• R(V) distribution (Gaussian)
• f_A distribution (Gaussian)
• star formation history (bins in log age)
• age-mass metallicity (single metallicity per log age?)

dust parameters should be coupled as our current physical model is for two sheets of dust - one in the foreground and one in the galaxy midplane. The A(V), R(V), and f_A distributions should be coupled such that the stars only seeing the foreground cloud should all be modeled by the 1st cloud A(V), R(V), and f_A distributions. The same for the stars behind the midplane cloud as they should see the combination of both clouds.

Need to have simulations to test and train the megabeast.
Start with A(V) only simulations that simulate the output of the BEAST. Later add in simulations of the BEAST inputs and run the full BEAST+MegaBEAST.

I've been using the make_naive_maps script for a while now, and I've always felt like there was something off about the coordinate system.

I believe I have now found the error. There's two mistakes when the RA and Dec of the sources are converted to pixel coordinates of the map WCS.

• The code loops over the pixels using range(n_x) and range(n_y), which start at 0. But inregions_for_objects, we have pixcrd = wcs_info.wcs_world2pix(world, 1), which asks for a 1-based index instead. This type of pixel coordinate is what's shown in DS9 for example. But when accessing the data array using numpy, 0-based coordinates are needed.
• The pixel coordinates in a fits image are relative to the center of the pixels, not the bottom left corner. The bottom left corner of pixel (1,1) is at (0.5, 0.5), and the top left corner at (1.5, 1.5) in case of 1-based coordinates. For 0-based, this pixel becomes (0, 0), with corners (-0.5, -0.5) and (0.5, 0.5) respectively.

So 1.6 would belong belongs to pixel 2, not pixel 1. Therefore, the floating point pixel coordinates should be rounded, not truncated.

# this is wrong
x = pixcrd[:, 0].astype(int)
# this should be ok
x = np.rint(pixcrd[:, 0]).astype(int)

This is fixed in my own branch, but there are some other changes. I'll put in a clean pull request once I'm more sure that my own code is correct.

Need to deal with this when reading them.
Pad the value array with zeros and give the index array any valid index (copy from early part). Zeros will mean that these padded indexes will not contribute to the integrated prob.

Inspired by @kmcquinn's ISM*-at-ST talk, the mixture model should have 3 terms. Stars internal to the galaxy being studied, the background galaxies/QSOs, and foreground stars.

The foreground stars can be modeled probabistically from a model of the Milky Way for the specific direction of the galaxy under study. This is exactly what is needed for the ensemble model.

The background galaxies could be determined from deep field observations (should include all usual filters).

By splitting the non-galaxy-being-studied model between foreground and background, it maybe more easy to generate and a higher fidelity model of the actual distribution of star fluxes.

The megabeast use cases will include fitting:

1. a single set of BEAST results with a single megabeast ensemble model
2. BEAST results organized into many pixels with a different megabeast ensemble model for each pixel

The megabeast core should support both cases. The 1st case will be especially useful for testing the megabest via simulations. It will also support ensemble fitting for clusters, etc.

(This is minor). Although the standard BEAST outputs are stored in a stats.fits file, referring to stats.fits in this line of make_naive_maps.py may result in an error if one has changed the name of the output file where stats.fits does not exist in a file name. (Actually the error I encountered ends up overwriting the very stats.fits file! )
For this reason I'd lobby for only referring to .fits instead of to stats.fits, unless values_per_pixel.hd5 without a preceding stats is an expected input somewhere else at this point and it's hard to change that.

Very useful to be able to limit the ensemble fitting to certain region of HR space (e.g, RGB).

Hello from Astropy Project!

The following updates to your package might be necessary to ensure compatibility with the latest stable version of Astropy:

• MPLBACKEND is now set to Agg in ci-helpers, packages expecting interactive plotting should override it in .travis.yml
• Astropy 3.1 is not compatible with Numpy <1.13 versions. If you want to keep testing these older Numpy versions, please use ASTROPY_VERSION=3.0 or ASTROPY_VERSION=LTS in your Travis CI matrix.
• Add sphinx-astropy as a package dependency if you are using astropy-helpers 3.1 or later. Otherwise, your documentation build (e.g., on ReadTheDocs) might fail.
• If you are using six that is bundled with Astropy, please consider using the standalone six package instead.

If these are no longer applicable to you, please close the issue.

This is an automated issue for packages that opted in for automated astropy-helpers update. If this is opened in error, please let @pllim know.

xref astropy/astropy-tools#108

Need to articulate the full ensemble model for the megaBEAST.
Dust properties.
Stellar properties.
Distance.

As I'm working on naive A_V maps (and this is likely relevant to MegaBEAST A_V maps too), I'm wondering if we've thought about the impact of the high A_V AGB failure mode. Specifically, what cuts, if any, should we make to avoid including them in the maps? If I recall correctly, they have good fit quality, so they wouldn't be eliminated with a chi^2 cut. When I was finding sources for my HST IC1613 spectroscopic UV extinction curve proposal, I had to remove BEAST fits with A_V >~ 2, but that was of course for a different purpose.

To whom it may concern,

If you are using https://github.com/astropy/ci-helpers in your appveyor.yml , please know that the Astropy project has dropped active development/support for Appveyor CI. If it still works, good for you, because we did not remove the relevant files (yet). But if it ever stops working, we have no plans to fix anything for Appveyor CI. Please consider using native Windows support other CI, e.g., Travis CI (see https://docs.travis-ci.com/user/reference/windows/). We apologize for any inconvenience caused.

If this issue is opened in error or irrelevant to you, feel free to close. Thank you.

xref astropy/ci-helpers#464

For the BEAST, it's straightforward to write a test that fits a single SED, but the MegaBEAST requires a bunch of fits. Is there PHAT data we could use for this? (I don't want to use METAL fits until they're published.) Or perhaps make some fake data? We're going to need to do that anyway for verification.