# Galaxy Modeling Plans for Summer 2013

## High-level Overview (read this first !!)

See the document attached here.

## Models & Free Parameters

• lensfit uses "bulge+disk" models, but with ellipticity and radius ratio fixed between components; parameters = 1 ellipse (3 or 5 parameters, depending on whether to fit centroid) + 2 amplitudes.
• All Sersic components will be approximated by sums of Gaussians, and convovled with multi-shapelet PSF. Sersic index for each component will be fixed.
• Will have support for fitting arbitrary number of linear components (two = infinity), but all components must have the same ellipticity and fixed radius ratio.
• Because it's all sums of Gaussians (with some amplitude ratios fixed and others free), can use any combination of Sersic components, or just fit single Gaussians.
• Current MultiShapelet modeling code does not include centroid as a free parameter, but we'll need to do that here (at least in the last iteration or two).

## Optimization & Sampling

• To simulate lensfit, we need to use Monte Carlo or grid sampling, not a greedy optimizer. Greedy optimizers are a pain anyway because they're much more mathematically complicated, and I don't fully trust the one I'm using.
• We might still want to use a greedy optimizer to seed the Monte Carlo...but that would also required getting covariance matrices out of the greedy optimizer, which I'm not doing at present, and could involve a significant amount of development (in that it might require adding derivatives to models).
• I'd like to start by trying to jump straight from moments (and errors on the moments) to adaptive Monte Carlo, as this would keep me from having to deal with optimizers or shapelet-model derivatives at all in S13.
• If we were to start with a greedy optimizer, it would only affect the performance in StackFit mode, and hence trying to skip it would not affect our measurement of the marginal cost of MultiFit.
• While it's quite possible we'll want to use a greedy optimizer to start with in the end, we may only gain from it if we have time to get it working well (time I don't think we have in S13).
• This reduces the overlap between the existing MultiShapelet model-fitting code and the work to be done in S13, which can be good (didn't trust that optimizer anyway) and bad (less synergy with HSC work on model mags).
• Regardless of whether I try to use a greedy optimizer at all, I will add an option to do a brute-force grid sampling, as this will be very important for debugging and testing the Monte Carlo sampler.

## Implementation

(class APIs t.b.d.)

• Though I initially planned to write this as another generic source measurement plugin (at least for the single-epoch fitter), I am now planning to write it as a separate task that would run after processCoadd.py and generate its own source table. This has a number of advantages:
• Easier to rerun just the code I am interested in when debugging.
• Will be able to implement the actual fitting in Python (with model evaluation still done in C++), which should greatly aid in debugging and visualization.
• If things go extremely poorly on the optimizer front, would be able to plug in SciPy optimizers easily.
• So far, model evaluation has always dominated the computation time, so we shouldn't take a significant performance hit. If we do, I can easily code it up again in C++ after debugging and tuning in Python.
• Multi-epoch fitter would need to be a new task regardless; easiest to hard-code this to the model fitting I'm doing rather than make a generic MultiFit driver task.
• By working more in Python (less in C++) and using Task framework, I could use existing mechanisms to do timing, but writing timing directly to the output catalog still has some advantages:
• Timing outputs will be large compared to task timings now, as there'll be many timings per object, and that may not work well with existing mechanisms.
• It'd be very convenient to have timing results in the same catalog as the fitted parameters, as that would let us easily plot them against each other.

## Work Plan

1. Update the model evaluation code to handle a combination of fixed and free amplitudes in the models.
• Instead of outputting an nPixels vector, needs to output an nPixels x nAmplitudes matrix.
• Need to use a matrix as input to set fixed amplitude ratios between Gaussian components and which linear combinations are to be free amplitudes.
1. Implement brute-force grid sampler and visualization of likelihood surfaces, in first version of single-epoch fitting Task.
• This may inform how to implement the adaptive Monte-Carlo sampler, and whether I need to reconsider trying to avoid starting with a greedy optimizer.
• Will need at least some simulations with representative S/N, resolution, and morphological distribution at this stage.
1. Add priors to model evaluation.
• Probably just use priors from lensfit papers for now.
• Real-world priors may not match what went into ImSim catalogs, so maybe we should get priors from that, and add scatter?
1. Implement adaptive Monte Carlo sampler in single-epoch mode.
1. Implement multi-epoch mode.