Enable multiband runs (single-band unchanged)

[Jayashree-Behera]

This PR adds multiband source fitting utilities to crowdsource_base.py and updates wise_proc.py to run in multiband mode. Multiband functions are named after their corresponding single-band functions with "_multiband" appended.
Single-band run behavior is unchanged.

[schlafly]

This looks great, thanks Shree.

Let's delete mosaic.py which will be broken. I think we will need some decam_proc.py updates but they should be modest.

I'll have some more comments tomorrow, sorry to be slow.

[schlafly]

You're duplicating my existing code, but I can't figure out for the life of me why I did this. Let's keep this for this commit, but can you file an issue to check removing this after we merge this commit?

[schlafly]

This looks good. Several comments inline. One bigger question: when running on a multiband image, I think fit_im now returns a "multiband" output; i.e., both _joint quantities and the original quantities, now with _b0 appended. There's also one newer quantity, snr_joint.

I'm wondering whether in the single band case we should compress the output before returning, doing roughly:
cat = {k[:-3]: cat[k] for k in cat if k.endswith('_b0')}
Then I think e.g. DECam processing and single-band WISE processing would be less affected (the returned catalogs would be basically the same as before. You may have already addressed the single-band WISE case in wise_proc?

[schlafly]

Looks like dead code.

[schlafly]

I think we lost a bit of the logic here. It clearly doesn't make a big impact, but here's the intent:

1. we cut the image into subregions
2. each subregion has a 'primary' region; the union of the primary regions covers the image.
3. it also has a boundary region; we include the boundary regions in the fits so that stars on the edge of the primary region aren't cut off and can be fit simultaneously with their neighbors. The primary + boundary region gets the 'a' suffix ('all')
4. we do the fit on each subregion independently ('all')
5. we only fill out the fluxes of the stars that were in the 'primary' region

I think the new code does all of this right except for 5.

[schlafly]

I think I personally would have done flux_avg = np.sum(flux * fluxunc**-2) / np.sum(fluxunc ** -2), and xcen /= flux_avg, ycen /= flux_avg. This is also fine though since it's the first iteration and we just need to do something sane.

Perhaps add some comments that this only applies to the first iteration since that's the iteration where we don't know the flux scaling.

What happens to the centroids for newly detected sources? It feels like they probably need similar handling?

[Jayashree-Behera]

The derivative columns in build_sparse_matrix are flux-scaled using guessflux from previous iteration, so the LSQR dx/dy coefficients already behave like centroid shifts. Except for the first iteration and so we do xcen /= flux_snr normalization.

[schlafly]

Suggested change

	guessmbda = np.concatenate([guess_local_flux, skypar.get((bdxf, bdyf), np.zeros(B*nskypar, 'f4'))]) if nskypar>0 else guess_local_flux
	guessmbda = np.concatenate([guess_local_flux, skypar[(bdxf, bdyf)]])

I feel like on 1164 you make sure that this exists so you don't need to guard against it not existing here?

[Jayashree-Behera]

Right, no need for both lines to do this. I will correct that.

[schlafly]

This got a little more involved than I immediately thought necessary. Would you mind looking at the following and seeing what I'm missing?

            if numpy.all(numpy.isnan(tmodel)):
                raise ValueError("Model is all NaNs")
            ind_pri_in_full = numpy.flatnonzero(mbd)
            ind_pri_in_all = numpy.flatnonzero(mbd[mbda])
            Nall = np.sum(mbda)

            # copy the centroid parameters for primary sources from the subregion fit
            # into the full image parameter set
            if tpsfderiv:
                for i in range(2):
                flux[N*B + ind_pri_in_full + i] = tflux[0][Nall*B + ind_pri_in_all + i]
                 
            skypar[(bdxf, bdyf)] = tflux[0][-B*nskypar:].copy()

            # copy the flux parameters for primary sources from the subregion fit
            # into the full image parameter set
            for b in range(B):
                models[b][spri] = tmodel[b][sfit]
                msky[b][spri] = tmsky[b][sfit]
                flux[ind_pri_in_full*B + b] = tflux[0][ind_pri_in_all*B + b]
                for k in range(n_psfcomps):
                    if Nall == 0:
                        continue
                    psf_stamps[b][k][ind_pri_in_full] = [
			 psfmod.central_stamp(psfsbda[i][tind], minsz)
                                    for tind in ind_pri_in_all]

(roughly I'm trying to follow the old logic more closely, but I'm probably missing something you're doing there!)

[Jayashree-Behera]

Thanks for writing it out. Yes, I originally tried to follow the old logic more closely. The reason I introduced the explicit ind_all to ind_pri_loc mapping is that the subregion fit is indexed in the ALL-sources ordering (flatnonzero(mbda_src)), and I wanted to be explicit about mapping PRIMARY sources back into that local ordering.

If we are confident that mbd[mbda] preserves exactly the same ordering as flatnonzero(mbda_src), then your version should indeed work and would be cleaner.

[Jayashree-Behera]

That checks out, so I made that change.

Regarding other parts of the block, the logic is much clearer in your form. It would definitely be nicer if the parameter layout were source-major like that, since the indexing becomes much shorter and cleaner.

But in my implementation, the flux vector is band-major (b*N + i) where b = band and i = source, so I had to index it slightly differently to stay consistent with how the matrix is built. If it’s not a strong preference, I’d lean toward keeping the current layout for now to avoid changing any of my analysis notebooks immediately.