Coverage for metadetect / fitting.py: 4%
630 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-09 01:50 -0800
1import logging
2import copy
4import numpy as np
6import ngmix
7import ngmix.prepsfmom
8from ngmix.gexceptions import (
9 BootPSFFailure, PSFFluxFailure,
10)
11from ngmix.moments import make_mom_result, fwhm_to_T
12from ngmix.bootstrap import bootstrap
13from ngmix.runners import Runner, PSFRunner
14from ngmix.guessers import SimplePSFGuesser
15from ngmix.fitting import Fitter
16from packaging.version import parse
18from .util import Namer
19from . import procflags
21MAX_NUM_SHEAR_BANDS = 6
23logger = logging.getLogger(__name__)
25if parse(ngmix.__version__) < parse("2.1.0"): 25 ↛ 26line 25 didn't jump to line 26 because the condition on line 25 was never true
26 MOMNAME = "mom"
27else:
28 MOMNAME = "sums"
31def fit_mbobs_gauss(
32 *,
33 mbobs,
34 bmask_flags,
35 rng,
36 shear_bands=None,
37 obj_runner=None,
38 psf_runner=None,
39 coadd=False,
40):
41 """Fit a multiband obs using a Gaussian fit.
43 Parameters
44 ----------
45 mbobs : ngmix.MultiBandObsList
46 The observation to use for shear measurement.
47 bmask_flags : int
48 Observations with these bits set in the bmask are not fit.
49 rng : np.random.RandomState
50 Random state for fitting.
51 shear_bands : list of int, optional
52 A list of indices into each mbobs that denotes which band is used for shear.
53 Default is to use all bands.
54 obj_runner : ngmix.runners.Runner, optional
55 If not None, the result of `get_gauss_obj_runner` is suggested. If None,
56 `get_gauss_obj_runner` is called.
57 psf_runner : ngmix.runners.PSFRunner, optional
58 If not None, the result of `get_gauss_psf_runner` is suggested. If None,
59 `get_gauss_psf_runner` is called.
60 coadd : bool, optional
61 If True, coadd the mbobs over all bands and then fit. Default is False.
63 Returns
64 -------
65 res : np.ndarray
66 A structured array of the fitting results.
67 """
68 if shear_bands is None:
69 shear_bands = list(range(len(mbobs)))
71 res = get_wavg_output_struct(len(mbobs), "gauss", shear_bands=shear_bands)
73 flags = 0
74 for obslist in mbobs:
75 if len(obslist) == 0:
76 flags |= procflags.MISSING_BAND
77 continue
79 if len(obslist) > 1:
80 flags |= procflags.INCONSISTENT_BANDS
81 continue
83 for obs in obslist:
84 if not np.any(obs.weight > 0):
85 flags |= procflags.ZERO_WEIGHTS
87 if np.any((obs.bmask & bmask_flags) != 0):
88 flags |= procflags.EDGE_HIT
90 if any(s >= len(mbobs) for s in shear_bands):
91 flags |= procflags.INCONSISTENT_BANDS
93 if coadd:
94 if flags == 0:
95 # first we coadd the shear bands
96 coadd_obs, coadd_flags = make_coadd_obs(mbobs, shear_bands=shear_bands)
97 flags |= coadd_flags
99 if flags == 0:
100 shear_mbobs = ngmix.observation.get_mb_obs(coadd_obs)
101 else:
102 if flags == 0:
103 shear_mbobs = ngmix.MultiBandObsList()
104 for band in shear_bands:
105 shear_mbobs.append(mbobs[band])
107 if flags == 0:
108 try:
109 ores = bootstrap(
110 shear_mbobs,
111 get_gauss_obj_runner(
112 rng,
113 len(shear_mbobs),
114 shear_mbobs[0][0].jacobian.get_scale(),
115 )
116 if obj_runner is None else obj_runner,
117 psf_runner=(
118 get_gauss_psf_runner(rng)
119 if psf_runner is None
120 else psf_runner
121 ),
122 )
123 except BootPSFFailure:
124 flags |= procflags.PSF_FAILURE
125 except PSFFluxFailure:
126 flags |= procflags.PSF_FAILURE
127 except Exception:
128 flags |= procflags.OBJ_FAILURE
129 flags |= procflags.PSF_FAILURE
131 if flags == 0:
132 res["gauss_obj_flags"] = ores["flags"]
133 res["gauss_T_flags"] = ores["flags"]
135 if len(res["gauss_band_flux_flags"].shape) == 2:
136 res["gauss_band_flux_flags"][0, shear_bands] = ores["flags"]
137 else:
138 res["gauss_band_flux_flags"] = ores["flags"]
140 if ores["flags"] == 0:
141 res["gauss_s2n"] = ores["s2n"]
142 res["gauss_g"] = ores["g"]
143 res["gauss_g_cov"] = ores["g_cov"]
144 res["gauss_T"] = ores["T"]
145 res["gauss_T_err"] = ores["T_err"]
147 if len(res["gauss_band_flux_flags"].shape) == 2:
148 res["gauss_band_flux"][0, shear_bands] = ores["flux"]
149 res["gauss_band_flux_err"][0, shear_bands] = ores["flux_err"]
150 else:
151 res["gauss_band_flux"] = ores["flux"]
152 res["gauss_band_flux_err"] = ores["flux_err"]
154 if coadd:
155 pres = shear_mbobs[0][0].psf.meta["result"]
156 res["gauss_psf_flags"] = pres["flags"]
157 if res["gauss_psf_flags"] == 0:
158 res["gauss_psf_T"] = pres["T"]
159 res["gauss_psf_g"] = _extract_g(pres)
160 if ores["flags"] == 0:
161 res["gauss_T_ratio"] = res["gauss_T"] / res["gauss_psf_T"]
162 else:
163 pflags = 0
164 psf_g_sum = np.zeros(2)
165 psf_T_sum = 0.0
166 wgt_sum = 0.0
167 for obslist in shear_mbobs:
168 for obs in obslist:
169 pflags |= obs.psf.meta["result"]["flags"]
170 if obs.psf.meta["result"]["flags"] == 0:
171 msk = obs.weight > 0
172 if not np.any(msk):
173 pflags |= procflags.ZERO_WEIGHTS
174 else:
175 _wgt = np.median(obs.weight[msk])
176 psf_T_sum += obs.psf.meta["result"]["T"] * _wgt
177 g = _extract_g(obs.psf.meta["result"])
178 psf_g_sum += (
179 g
180 * _wgt
181 * obs.psf.meta["result"]["T"]
182 )
183 wgt_sum += _wgt
185 res["gauss_psf_flags"] = pflags
186 if res["gauss_psf_flags"] == 0:
187 res["gauss_psf_T"] = psf_T_sum / wgt_sum
188 res["gauss_psf_g"] = psf_g_sum / psf_T_sum
189 if ores["flags"] == 0:
190 res["gauss_T_ratio"] = res["gauss_T"] / res["gauss_psf_T"]
192 res["gauss_flags"] = res["gauss_obj_flags"] | res["gauss_psf_flags"]
193 else:
194 res["gauss_flags"] |= flags
196 return res
199def _extract_g(res):
200 """
201 We may have fit the PSF with adaptive moments, so check for
202 'e' as well
203 """
204 if 'g' in res:
205 g = res['g']
206 elif 'e' in res:
207 e1, e2 = res['e']
208 g1, g2 = ngmix.shape.e1e2_to_g1g2(e1, e2)
209 g = np.array([g1, g2])
210 return g
213def get_gauss_psf_runner(rng):
214 psf_guesser = SimplePSFGuesser(
215 rng=rng,
216 guess_from_moms=True,
217 )
218 psf_fitter = Fitter(
219 model="gauss",
220 fit_pars={
221 "maxfev": 2000,
222 "xtol": 1.0e-5,
223 "ftol": 1.0e-5,
224 },
225 )
226 psf_runner = PSFRunner(
227 fitter=psf_fitter,
228 guesser=psf_guesser,
229 ntry=2,
230 )
231 return psf_runner
234def get_gauss_obj_runner(rng, nband, scale):
235 prior = _make_ml_prior(rng, scale, nband)
237 fitter = ngmix.fitting.Fitter(
238 model='gauss',
239 prior=prior,
240 fit_pars={
241 "maxfev": 2000,
242 "xtol": 5.0e-5,
243 "ftol": 5.0e-5,
244 },
245 )
246 guesser = ngmix.guessers.TPSFFluxGuesser(
247 rng=rng,
248 T=0.25,
249 prior=prior,
250 )
251 runner = Runner(
252 fitter=fitter,
253 guesser=guesser,
254 ntry=2,
255 )
257 return runner
260def _make_ml_prior(rng, scale, nband):
261 """make the prior for the fitter.
263 Parameters
264 ----------
265 rng: np.random.RandomState
266 The random number generator
267 scale: float
268 Pixel scale
269 nband: int
270 number of bands
271 """
272 g_prior = ngmix.priors.GPriorBA(sigma=0.3, rng=rng)
273 cen_prior = ngmix.priors.CenPrior(
274 cen1=0, cen2=0, sigma1=scale, sigma2=scale, rng=rng,
275 )
276 T_prior = ngmix.priors.TwoSidedErf(
277 minval=-10.0,
278 width_at_min=0.03,
279 maxval=1.0e6,
280 width_at_max=1.0e5,
281 rng=rng,
282 )
283 F_prior = ngmix.priors.TwoSidedErf(
284 minval=-1.0e4,
285 width_at_min=1.0,
286 maxval=1.0e9,
287 width_at_max=0.25e8,
288 rng=rng,
289 )
290 F_prior = [F_prior] * nband
292 prior = ngmix.joint_prior.PriorSimpleSep(
293 cen_prior=cen_prior,
294 g_prior=g_prior,
295 T_prior=T_prior,
296 F_prior=F_prior,
297 )
299 return prior
302def fit_mbobs_list_joint(
303 *, mbobs_list, fitter_name, bmask_flags, rng, shear_bands=None,
304 symmetrize=True, coadd=False,
305):
306 """Fit the ojects in a list of ngmix.MultiBandObsList using a joint fitter.
308 The fitter is run per object on the bands used for shear.
310 Parameters
311 ----------
312 mbobs_list : a list of ngmix.MultiBandObsList
313 The observations to use for shear measurement.
314 fitter_name : str
315 The name of the fitter to use.
316 bmask_flags : int
317 Observations with these bits set in the bmask are not fit.
318 shear_bands : list of int, optional
319 A list of indices into each mbobs that denotes which band is used for shear.
320 Default is to use all bands.
321 rng : np.random.RandomState
322 Random state for fitting.
323 symmetrize : bool, optional
324 If True, apply 4-fold symmetry to the mask+weight map. Default is True.
325 Ignored for gaussian fits.
326 coadd : bool, optional
327 If True, coadd the mbobs over all bands and then fit. Default is False.
328 Ignored for adaptive moments which always coadds.
330 Returns
331 -------
332 res : np.ndarray
333 A structured array of the fitting results.
334 """
335 if fitter_name in ["am", "admom"]:
336 fit_func = fit_mbobs_admom
337 kwargs = {"runner": get_admom_runner(rng), 'symmetrize': symmetrize}
338 elif fitter_name == "gauss":
339 fit_func = fit_mbobs_gauss
340 kwargs = {"coadd": coadd}
341 if len(mbobs_list) > 0 and len(mbobs_list[0]) > 0 and len(mbobs_list[0][0]) > 0:
342 scale = mbobs_list[0][0][0].jacobian.get_scale()
343 if coadd:
344 nband = 1
345 else:
346 if shear_bands is None:
347 nband = len(mbobs_list[0])
348 else:
349 nband = len(shear_bands)
350 kwargs["obj_runner"] = get_gauss_obj_runner(rng, nband, scale)
351 kwargs["psf_runner"] = get_gauss_psf_runner(rng)
352 else:
353 raise RuntimeError("Joint fitter '%s' not recognized!" % fitter_name)
355 res = []
356 for i, mbobs in enumerate(mbobs_list):
357 _res = fit_func(
358 mbobs=mbobs,
359 bmask_flags=bmask_flags,
360 shear_bands=shear_bands,
361 rng=rng,
362 **kwargs,
363 )
364 res.append(_res)
366 if len(res) > 0:
367 return np.hstack(res)
368 else:
369 return None
372def get_admom_runner(rng):
373 fitter = ngmix.admom.AdmomFitter(rng=rng)
374 guesser = ngmix.guessers.GMixPSFGuesser(
375 rng=rng, ngauss=1, guess_from_moms=True,
376 )
377 return ngmix.runners.Runner(
378 fitter=fitter, guesser=guesser,
379 ntry=2,
380 )
383def fit_mbobs_admom(
384 *,
385 mbobs,
386 bmask_flags,
387 rng,
388 shear_bands=None,
389 runner=None,
390 symmetrize=True,
391):
392 """Fit a multiband obs using adaptive moments.
394 This function forms a coadd of the shear bands and then runs
395 adaptive moments on the coadd.
397 Parameters
398 ----------
399 mbobs : ngmix.MultiBandObsList
400 The observation to use for shear measurement.
401 bmask_flags : int
402 Observations with these bits set in the bmask are not fit.
403 rng : np.random.RandomState
404 Random state for fitting.
405 shear_bands : list of int, optional
406 A list of indices into each mbobs that denotes which band is used for shear.
407 Default is to use all bands.
408 runner : ngmix.runners.Runner, optional
409 If not None, the result of `get_admom_runner` is suggested. If None,
410 `get_admom_runner` is called.
411 symmetrize : bool, optional
412 If True, apply 4-fold symmetry to the mask+weight map. Default is True.
414 Returns
415 -------
416 res : np.ndarray
417 A structured array of the fitting results.
418 """
419 runner = get_admom_runner(rng) if runner is None else runner
420 nband = len(mbobs)
421 if shear_bands is None:
422 shear_bands = list(range(len(mbobs)))
423 res = get_wavg_output_struct(nband, "am", shear_bands=shear_bands)
425 flags = 0
426 for obslist in mbobs:
427 if len(obslist) == 0:
428 flags |= procflags.MISSING_BAND
429 continue
431 if len(obslist) > 1:
432 flags |= procflags.INCONSISTENT_BANDS
433 continue
435 for obs in obslist:
436 if not np.any(obs.weight > 0):
437 flags |= procflags.ZERO_WEIGHTS
439 if np.any((obs.bmask & bmask_flags) != 0):
440 flags |= procflags.EDGE_HIT
442 if flags == 0:
443 # first we coadd the shear bands
444 coadd_obs, coadd_flags = make_coadd_obs(mbobs, shear_bands=shear_bands)
445 flags |= coadd_flags
447 if flags == 0:
448 # then fit the PSF
449 try:
450 pres = runner.go(coadd_obs.psf)
451 except Exception:
452 flags |= procflags.PSF_FAILURE
453 else:
454 res["am_psf_flags"] = pres["flags"]
455 if pres["flags"] == 0:
456 res["am_psf_g"] = pres["e"]
457 res["am_psf_T"] = pres["T"]
459 if flags == 0:
460 # then fit the object
461 if symmetrize:
462 sym_coadd_obs = symmetrize_obs_weights(coadd_obs)
463 else:
464 sym_coadd_obs = coadd_obs
465 try:
466 gres = runner.go(sym_coadd_obs)
467 except Exception:
468 flags |= procflags.OBJ_FAILURE
469 # replace no attempt
470 res["am_flags"] = flags
471 else:
472 res["am_T_flags"] = gres["T_flags"]
473 if gres["T_flags"] == 0:
474 res["am_T"] = gres["T"]
475 res["am_T_err"] = gres["T_err"]
476 if pres["flags"] == 0:
477 res["am_T_ratio"] = res["am_T"] / res["am_psf_T"]
479 res["am_obj_flags"] = gres["flags"]
480 if gres["flags"] == 0:
481 res["am_s2n"] = gres["s2n"]
482 res["am_g"] = gres["e"]
483 res["am_g_cov"] = gres["e_cov"]
485 # this replaces the flags so they are zero and unsets the default of
486 # no attempt
487 res["am_flags"] = (res["am_psf_flags"] | res["am_obj_flags"])
489 else:
490 # this branch ensures noattempt remains set
491 res["am_flags"] |= flags
493 return res
496def make_coadd_obs(mbobs, shear_bands=None):
497 """Coadd the observations in an ngmix mbobs assuming they all have the same
498 shaped images and same Jacobians.
500 Parameters
501 ----------
502 mbobs : ngmix.MultibandObsList
503 The observation to use for shear measurement.
504 shear_bands : list of int, optional
505 A list of indices into each mbobs that denotes which band is used for shear.
506 Default is to use all bands.
508 Returns
509 -------
510 coadd_obs : ngmix.Observation
511 The coadded observation with the PSF set.
512 flags : int
513 Any flags for errors in making the coadd obs.
514 """
515 nbands = len(mbobs)
516 flags = 0
517 if shear_bands is None:
518 shear_bands = list(range(nbands))
520 # if we have one input band, then simply
521 # return if if the data is consistent
522 if nbands == 1:
523 if shear_bands[0] == 0:
524 return mbobs[0][0], flags
525 else:
526 flags |= procflags.INCONSISTENT_BANDS
527 return None, flags
529 # if we ask for a single band for shear, then return that band
530 # if the data is consistent
531 if len(shear_bands) == 1:
532 if shear_bands[0] < nbands:
533 return mbobs[shear_bands[0]][0], flags
534 else:
535 flags |= procflags.INCONSISTENT_BANDS
536 return None, flags
538 # at this point, we are coadding more than one band
539 # let's check for consistent data
541 # we can't use shear bands we do not have
542 if any(s >= nbands for s in shear_bands):
543 flags |= procflags.INCONSISTENT_BANDS
544 return None, flags
546 # all of the image shapes and jacobians have to be the same
547 fobs = mbobs[shear_bands[0]][0]
548 for shear_band in shear_bands[1:]:
549 obs = mbobs[shear_band][0]
550 if (
551 (repr(fobs.jacobian) != repr(obs.jacobian))
552 or (fobs.image.shape != obs.image.shape)
553 or (not fobs.has_psf())
554 or (not obs.has_psf())
555 or (repr(fobs.psf.jacobian) != repr(obs.psf.jacobian))
556 or (fobs.psf.image.shape != obs.psf.image.shape)
557 ):
558 flags |= procflags.INCONSISTENT_BANDS
559 return None, flags
561 # if we get here, we need to coadd
562 # first make the weights
563 wgts = np.zeros(len(shear_bands))
564 for i, band in enumerate(shear_bands):
565 obs = mbobs[band][0]
566 msk = obs.weight > 0
567 if not np.any(msk):
568 flags |= procflags.ZERO_WEIGHTS
569 # let's not waste time since we are missing a band
570 return None, flags
571 else:
572 wgts[i] = np.median(obs.weight[msk])
573 wgts /= np.sum(wgts)
575 # finish coadding
576 image = np.zeros_like(mbobs[shear_bands[0]][0].image)
577 weight = np.zeros_like(mbobs[shear_bands[0]][0].weight)
578 mfrac = np.zeros_like(mbobs[shear_bands[0]][0].image)
579 nmfrac = 0
580 noise = np.zeros_like(mbobs[shear_bands[0]][0].image)
581 nnoise = 0
582 bmask = np.zeros_like(mbobs[shear_bands[0]][0].image, dtype=np.int32)
583 nbmask = 0
584 ormask = np.zeros_like(mbobs[shear_bands[0]][0].image, dtype=np.int32)
585 normask = 0
586 meta = {}
587 psf_meta = {}
588 psf_image = np.zeros_like(mbobs[shear_bands[0]][0].psf.image)
590 for i, shear_band in enumerate(shear_bands):
591 obs = mbobs[shear_band][0]
592 wgt = wgts[i]
594 meta.update(obs.meta)
595 meta.update(mbobs[shear_band].meta)
596 psf_meta.update(obs.psf.meta)
598 psf_image += (wgt * obs.psf.image)
600 image += (wgt * obs.image)
601 weight += (wgt**2 / obs.weight)
603 if obs.has_mfrac():
604 mfrac += (wgt * obs.mfrac)
605 nmfrac += 1
607 if obs.has_noise():
608 noise += (wgt * obs.noise)
609 nnoise += 1
611 if obs.has_bmask():
612 bmask |= obs.bmask
613 nbmask += 1
615 if obs.has_ormask():
616 ormask |= obs.ormask
617 normask += 1
619 weight = 1.0 / weight
620 msk = ~np.isfinite(weight)
621 if np.any(msk):
622 weight[msk] = 0
624 if np.all(weight == 0):
625 flags |= procflags.ZERO_WEIGHTS
626 return None, flags
628 for var in [nmfrac, nnoise, nbmask, normask]:
629 if var not in [0, len(shear_bands)]:
630 flags |= procflags.INCONSISTENT_BANDS
631 return None, flags
633 meta.update(mbobs.meta)
635 # set non-zero weights for PSF
636 psf_nse = np.sqrt(np.sum(psf_image**2)) / 1000.0
637 psf_weight = psf_image*0 + 1.0/psf_nse**2
639 kwargs = {
640 "weight": weight,
641 "jacobian": mbobs[shear_bands[0]][0].jacobian,
642 "psf": ngmix.Observation(
643 psf_image,
644 jacobian=mbobs[shear_bands[0]][0].psf.jacobian,
645 meta=psf_meta,
646 weight=psf_weight,
647 ),
648 "meta": meta,
649 }
650 if nmfrac > 0:
651 kwargs["mfrac"] = mfrac
652 if nnoise > 0:
653 kwargs["noise"] = noise
654 if nbmask > 0:
655 kwargs["bmask"] = bmask
656 if normask > 0:
657 kwargs["ormask"] = ormask
659 cobs = ngmix.Observation(image, **kwargs)
661 return cobs, flags
664def get_coellip_ngauss(name):
665 ngauss = int(name[7:])
666 return ngauss
669def fit_all_psfs(mbobs, rng):
670 """
671 measure all psfs in the input observations and store the results
672 in the meta dictionary, and possibly as a gmix for model fits
674 Parameters
675 ----------
676 mbobs: ngmix.MultiBandObsList
677 The observations to fit
678 psf_conf: dict
679 Config for the measurements/fitting
680 rng: np.random.RandomState
681 The random number generator, used for guessers
682 """
683 fitter = ngmix.admom.AdmomFitter(rng=rng)
684 guesser = ngmix.guessers.GMixPSFGuesser(
685 rng=rng, ngauss=1, guess_from_moms=True,
686 )
688 runner = ngmix.runners.PSFRunner(
689 fitter=fitter, guesser=guesser, ntry=10,
690 )
692 for obslist in mbobs:
693 assert len(obslist) == 1, 'metadetect is not multi-epoch'
695 obs = obslist[0]
696 runner.go(obs=obs)
698 flags = obs.psf.meta['result']['flags']
699 if flags != 0:
700 raise BootPSFFailure("failed to measure psfs: %s" % flags)
703def fit_mbobs_list_wavg(
704 *, mbobs_list, fitter, bmask_flags, shear_bands=None, fwhm_reg=0,
705 symmetrize=True,
706):
707 """Fit the ojects in a list of ngmix.MultiBandObsList using a weighted average
708 over bands.
710 The fitter is run per-band and then results are combined across bands
711 via a weighted average.
713 Parameters
714 ----------
715 mbobs_list : a list of ngmix.MultiBandObsList
716 The observations to use for shear measurement.
717 fitter : ngmix fitter
718 The fitter to use per band per MultiBandObsList.
719 bmask_flags : int
720 Observations with these bits set in the bmask are not fit.
721 shear_bands : list of int, optional
722 A list of indices into each mbobs that denotes which band is used for shear.
723 Default is to use all bands.
724 fwhm_reg : float, optional
725 This value is converted to T and used to regularize the shapes via
727 e_{1,2} = M_{1,2}/(T + T_reg)
729 For Gaussians, this relationship is equivalent to smoothing by a round
730 Gaussian with FWHM `fwhm_reg`.
731 symmetrize : bool, optional
732 If True, apply 4-fold symmetry to the mask+weight map. Default is True.
734 Returns
735 -------
736 res : np.ndarray
737 A structured array of the fitting results.
738 """
739 res = []
740 for i, mbobs in enumerate(mbobs_list):
742 _res = fit_mbobs_wavg(
743 mbobs=mbobs,
744 fitter=fitter,
745 bmask_flags=bmask_flags,
746 shear_bands=shear_bands,
747 fwhm_reg=fwhm_reg,
748 symmetrize=symmetrize,
749 )
750 res.append(_res)
752 if len(res) > 0:
753 return np.hstack(res)
754 else:
755 return None
758def fit_mbobs_wavg(
759 *,
760 mbobs,
761 fitter,
762 bmask_flags,
763 shear_bands=None,
764 fwhm_reg=0,
765 symmetrize=True,
766):
767 """Fit the object in the ngmix.MultiBandObsList using a weighted average
768 over bands.
770 The fitter is run per-band and then results are combined across bands
771 via a weighted average.
773 Parameters
774 ----------
775 mbobs : ngmix.MultiBandObsList
776 The observation to use for shear measurement.
777 fitter : ngmix fitter
778 The fitter to use per band.
779 bmask_flags : int
780 Observations with these bits set in the bmask are not fit.
781 shear_bands : list of int, optional
782 A list of indices into each mbobs that denotes which band is used for shear.
783 Default is to use all bands.
784 fwhm_reg : float, optional
785 This value is converted to T and used to regularize the shapes via
787 e_{1,2} = M_{1,2}/(T + T_reg)
789 For Gaussians, this relationship is equivalent to smoothing by a round
790 Gaussian with FWHM `fwhm_reg`.
791 symmetrize : bool, optional
792 If True, apply 4-fold symmetry to the mask+weight map. Default is True.
794 Returns
795 -------
796 res : np.ndarray
797 A structured array of the fitting results.
798 """
799 nband = len(mbobs)
800 all_res = []
801 all_psf_res = []
802 all_is_shear_band = []
803 all_wgts = []
804 all_flags = []
806 if shear_bands is None:
807 shear_bands = list(range(nband))
809 if fitter.kind == 'am':
810 assert len(mbobs) == 1, 'Use only one band for adaptive moments'
812 for band in range(nband):
813 all_is_shear_band.append(True if band in shear_bands else False)
814 fres = _fit_obslist(
815 obslist=mbobs[band],
816 fitter=fitter,
817 bmask_flags=bmask_flags,
818 symmetrize=symmetrize,
819 )
820 all_wgts.append(fres["wgt"])
821 all_res.append(fres["obj_res"])
822 all_psf_res.append(fres["psf_res"])
823 all_flags.append(fres["flags"])
825 return _combine_fit_results_wavg(
826 all_res=all_res,
827 all_psf_res=all_psf_res,
828 all_is_shear_band=all_is_shear_band,
829 all_wgts=all_wgts,
830 model=fitter.kind,
831 all_flags=all_flags,
832 shear_bands=shear_bands,
833 fwhm_reg=fwhm_reg,
834 )
837def _fit_obslist(
838 *,
839 obslist,
840 fitter,
841 bmask_flags,
842 symmetrize,
843):
844 if len(obslist) == 0:
845 # we will flag this later
846 res = {}
847 res["flags"] = 0
848 res["flags"] |= procflags.MISSING_BAND
849 res["wgt"] = 0
850 res["obj_res"] = None
851 res["psf_res"] = None
852 return res
853 else:
854 return _fit_obs(
855 obs=obslist[0],
856 fitter=fitter,
857 bmask_flags=bmask_flags,
858 symmetrize=symmetrize,
859 )
862def _fit_obs(
863 *,
864 obs,
865 fitter,
866 bmask_flags,
867 symmetrize,
868):
869 if isinstance(fitter, ngmix.prepsfmom.PrePSFMom):
870 psf_go_kwargs = {"no_psf": True}
871 else:
872 psf_go_kwargs = {}
874 res = {}
875 flags = 0
877 if symmetrize:
878 obs = symmetrize_obs_weights(obs)
880 if not np.any(obs.weight > 0):
881 flags |= procflags.ZERO_WEIGHTS
883 if np.any((obs.bmask & bmask_flags) != 0):
884 flags |= procflags.EDGE_HIT
886 if flags != 0:
887 # we will flag this later
888 res["flags"] = flags
889 res["wgt"] = 0
890 res["obj_res"] = None
891 res["psf_res"] = None
892 else:
893 res["flags"] = flags
894 res["wgt"] = np.median(obs.weight[obs.weight > 0])
895 res["obj_res"] = fitter.go(obs)
896 res["psf_res"] = fitter.go(obs.psf, **psf_go_kwargs)
898 if fitter.kind == "am" and MOMNAME == "mom":
899 res["obj_res"]["mom"] = res["obj_res"]["sums"]
900 res["obj_res"]["mom_cov"] = res["obj_res"]["sums_cov"]
901 res["psf_res"]["mom"] = res["psf_res"]["sums"]
902 res["psf_res"]["mom_cov"] = res["psf_res"]["sums_cov"]
904 if res["obj_res"]["flags"] != 0:
905 logger.debug("per band fitter failed: %s" % res["obj_res"]['flagstr'])
907 if res["psf_res"]["flags"] != 0:
908 logger.debug("per band psf fitter failed: %s" % res["psf_res"]['flagstr'])
910 return res
913def _sum_bands_wavg(
914 *, all_res, all_is_shear_band, all_wgts, all_flags, all_wgt_res,
915):
916 """A function to sum all of the moments across different bands for combined moments
917 estimates.
919 This function sums the moments together using the entires from `all_wgts`.
921 raw_mom = \\sum_i all_wgts[i] * all_res[i]["sums"]
923 If `all_wgt_res` is given, then the weights are adjust to include the ratio of the
924 flux moments between `all_wgt_res` and `all_res`:
926 all_wgts[i] -> all_wgts[i] * all_wgt_res[i]["sums"][5] / all_res[i]["sums"][5]
928 If `mom_norm` is in the results dicts in `all_res` and `all_wgt_res`, then the
929 moments fields `mom` are divided by this normalization before being averaged.
931 The fluxes are summed without any additional weighting beyond the entires in
932 `all_wgts`
934 flux = \\sum_i all_wgts[i] * all_res[i]["flux"]
937 Similar sums are carried out for the the moments covariance and the flux variance
938 with the approprtiate weights squared.
940 Parameters
941 ----------
942 all_res : list of dicts
943 List of the moments fit results from each band.
944 all_is_shear_band : list of bool
945 List of bools indicating if a given band is to be summed into the outputs and
946 thus used for shear.
947 all_wgts : list of floats
948 List of floats giving the total weight for each band.
949 all_flags : list of ints
950 List of ints indicating if any of the input bands are flagged for some reason.
951 These flags are ORed into the final_flags outputs.
952 all_wgt_res : list of dicts, optional
953 If not None, then additional weighting by the flux moment ratio of
954 `all_wgt_res[i]["sums"][5] / all_res[i]["sums"][5]` is applied to the weighting
955 for the moments sums. This field is used to weight the PSF moment sums by
956 the fluxes of the objects so that the band-averaged sizes of stars match the
957 band-averaged PSF size.
959 Returns
960 -------
961 res : dict
962 A dictionary with the following keys:
964 raw_mom : np.ndarray
965 The array of summed moments. You need to divide by `wgt_sum` to
966 normalize the moments properly.
967 raw_mom_cov : np.ndarray
968 The covariance of the summed moments. You need to divide by
969 `wgt_sum**2` to normalize the moments properly.
970 wgt_sum : float
971 The sum of the weights used in the moments sum. To get final moments,
972 divide `raw_mom` by `wgt_sum` and `raw_mom_cov` by `wgt_sum**2`.
973 final_flags : int
974 A final set of bit flags for the summed moments. This field marks
975 missing bands, missing moments, or zero weights. It is also ORed
976 with all of the input `all_flags` entires.
977 used_shear_bands : list of bool
978 A list of which bands were used for shear. This field can be used to
979 track if two successive calls to this function with slightly different
980 inputs end up using the same bands.
981 flux : float
982 The total flux. You must divide by `wgt_sum` to get the output
983 flux.
984 flux_var : float
985 The variance in the total flux. You must divide by `wgt_sum**2`
986 to get the output flux variance.
987 """
988 tot_nband = len(all_res)
989 raw_mom = np.zeros(6, dtype=np.float64)
990 raw_mom_cov = np.zeros((6, 6), dtype=np.float64)
991 wgt_sum = 0.0
992 used_shear_bands = [False] * tot_nband
993 final_flags = 0
994 flux = 0.0
995 flux_var = 0.0
997 for iband, (wgt, res, issb, flags) in enumerate(zip(
998 all_wgts, all_res, all_is_shear_band, all_flags
999 )):
1000 if all_wgt_res is not None:
1001 wgt_res = all_wgt_res[iband]
1002 else:
1003 wgt_res = None
1005 if issb:
1006 # the input flags mark very basic failures and are ORed across all bands
1007 # these are things like missing and or all zero-weight data, edges, etc.
1008 final_flags |= flags
1010 # we mark missing data or moments for PSF and objects separately
1011 if res is None or (all_wgt_res is not None and wgt_res is None):
1012 final_flags |= procflags.MISSING_BAND
1013 elif (
1014 (MOMNAME not in res or MOMNAME+"_cov" not in res)
1015 or (
1016 all_wgt_res is not None and wgt_res is not None and (
1017 MOMNAME not in wgt_res or MOMNAME+"_cov" not in wgt_res
1018 )
1019 )
1020 ):
1021 final_flags |= procflags.NOMOMENTS_FAILURE
1023 if (
1024 res is not None
1025 and MOMNAME+"_norm" in res
1026 and np.isfinite(res[MOMNAME+"_norm"])
1027 ):
1028 mom_norm = res[MOMNAME+"_norm"]
1029 else:
1030 mom_norm = 1.0
1032 if (
1033 all_wgt_res is not None
1034 and wgt_res is not None
1035 and res is not None
1036 and MOMNAME in res
1037 and MOMNAME in wgt_res
1038 ):
1039 if res[MOMNAME][5] != 0:
1040 if (
1041 MOMNAME+"_norm" in res
1042 and np.isfinite(res[MOMNAME+"_norm"])
1043 and MOMNAME+"_norm" in wgt_res
1044 and np.isfinite(wgt_res[MOMNAME+"_norm"])
1045 ):
1046 if wgt_res[MOMNAME+"_norm"] == 0:
1047 flux_mom_ratio = 1.0
1048 final_flags |= procflags.ZERO_WEIGHTS
1049 else:
1050 flux_mom_ratio = (
1051 wgt_res[MOMNAME][5]
1052 / wgt_res[MOMNAME+"_norm"]
1053 / res[MOMNAME][5]
1054 * res[MOMNAME+"_norm"]
1055 )
1056 else:
1057 flux_mom_ratio = wgt_res[MOMNAME][5] / res[MOMNAME][5]
1058 else:
1059 flux_mom_ratio = 1.0
1060 final_flags |= procflags.ZERO_WEIGHTS
1061 else:
1062 # we do not flag zero weight here since this is a missing band
1063 # or missign moments or no flux weighting
1064 flux_mom_ratio = 1.0
1066 if wgt <= 0:
1067 final_flags |= procflags.ZERO_WEIGHTS
1069 if res is not None and MOMNAME in res and MOMNAME+"_cov" in res:
1070 flux += (wgt * res[MOMNAME][5])
1071 flux_var += (wgt**2 * res[MOMNAME+"_cov"][5, 5])
1073 # there are a few factors here
1074 # wgt - the averaging weight for the moments
1075 # flux_mom_ratio - the ratio of the flux in wgt_res to res
1076 # This factor used to properly weight PSF model averages for objects
1077 # since the sums are actually things like flux * T. So for the PSF T
1078 # we actually average
1079 # flux_mom_ratio * sums[4] = (object flux/psf flux) * (psf_flux * T)
1080 # since sums[4] = flux * T.
1081 # mom_norm - sum of the moments weight function
1082 # This factor removes the moments dependence on area of the stamp.
1083 # The flux averages above do not get these factors since we want the
1084 # weight function to peak at 1 for flux measurements (and only care
1085 # about the object).
1087 # we added an extra sum for AM, for rho4
1088 rmoms = res[MOMNAME][:raw_mom.size]
1089 rcovs = res[MOMNAME + '_cov'][:raw_mom.size, :raw_mom.size]
1091 raw_mom += (wgt * flux_mom_ratio / mom_norm * rmoms)
1092 raw_mom_cov += (
1093 (wgt * flux_mom_ratio / mom_norm)**2
1094 * rcovs
1095 )
1096 wgt_sum += wgt
1098 used_shear_bands[iband] = True
1100 # make sure we flag missing data or all zero weight sums
1101 if sum(used_shear_bands) > 0 and wgt_sum <= 0:
1102 final_flags |= procflags.ZERO_WEIGHTS
1104 return dict(
1105 raw_mom=raw_mom,
1106 raw_mom_cov=raw_mom_cov,
1107 wgt_sum=wgt_sum,
1108 final_flags=final_flags,
1109 used_shear_bands=used_shear_bands,
1110 flux=flux,
1111 flux_var=flux_var,
1112 )
1115def _make_mom_res(*, raw_mom, raw_mom_cov, raw_flux, raw_flux_var, fwhm_reg):
1116 if fwhm_reg > 0:
1117 momres_t = make_mom_result(raw_mom, raw_mom_cov)
1119 T_reg = fwhm_to_T(fwhm_reg)
1121 # the moments are not normalized and are sums, so convert T_reg to a sum using
1122 # the flux sum first via T_reg -> T_reg * raw_mom[5]
1123 amat = np.eye(6)
1124 amat[4, 5] = T_reg
1126 raw_mom_orig = raw_mom.copy()
1127 if np.isnan(raw_mom_orig[0]):
1128 raw_mom[0] = 0
1129 if np.isnan(raw_mom_orig[1]):
1130 raw_mom[1] = 0
1131 reg_mom = np.dot(amat, raw_mom)
1132 if np.isnan(raw_mom_orig[0]):
1133 raw_mom[0] = np.nan
1134 reg_mom[0] = np.nan
1135 if np.isnan(raw_mom_orig[1]):
1136 raw_mom[1] = np.nan
1137 reg_mom[1] = np.nan
1139 reg_mom_cov = np.dot(amat, np.dot(raw_mom_cov, amat.T))
1140 momres = make_mom_result(reg_mom, reg_mom_cov)
1142 # use old T
1143 for col in ["T", "T_err", "T_flags", "T_flagstr"]:
1144 momres[col] = momres_t[col]
1145 momres["flags"] |= momres_t["flags"]
1146 else:
1147 momres = make_mom_result(raw_mom, raw_mom_cov)
1149 momres["flux"] = raw_flux
1150 if momres["flux"] <= 0:
1151 momres["flags"] |= ngmix.flags.NONPOS_FLUX
1153 if raw_flux_var > 0:
1154 momres["flux_err"] = np.sqrt(raw_flux_var)
1155 momres["s2n"] = momres["flux"] / momres["flux_err"]
1156 else:
1157 momres["flux_flags"] |= ngmix.flags.NONPOS_VAR
1158 momres["flux_err"] = np.nan
1159 momres["s2n"] = np.nan
1160 momres["flags"] |= ngmix.flags.NONPOS_VAR
1162 momres["flux_flagstr"] = procflags.get_procflags_str(momres["flux_flags"])
1163 momres["flagstr"] = procflags.get_procflags_str(momres["flags"])
1165 return momres
1168def _combine_fit_results_wavg(
1169 *, all_res, all_psf_res, all_is_shear_band, all_wgts, model, all_flags, shear_bands,
1170 fwhm_reg,
1171):
1172 tot_nband = len(all_res)
1173 nband = (
1174 sum(1 if issb else 0 for issb in all_is_shear_band)
1175 if all_is_shear_band
1176 else 0
1177 )
1178 nonshear_nband = (
1179 sum(0 if issb else 1 for issb in all_is_shear_band)
1180 if all_is_shear_band
1181 else 0
1182 )
1183 blens = [
1184 len(all_res), len(all_psf_res), len(all_is_shear_band),
1185 len(all_wgts), len(all_flags)
1186 ]
1188 n = Namer(front=model)
1190 data = get_wavg_output_struct(tot_nband, model, shear_bands=shear_bands)
1192 if (
1193 tot_nband == 0
1194 or nband == 0
1195 or tot_nband != nband + nonshear_nband
1196 or not all(b == tot_nband for b in blens)
1197 ):
1198 if nband == 0 or tot_nband == 0:
1199 psf_flags = procflags.MISSING_BAND
1200 mdet_flags = procflags.MISSING_BAND
1201 flux_flags = [procflags.MISSING_BAND] * tot_nband
1202 else:
1203 psf_flags = procflags.INCONSISTENT_BANDS
1204 mdet_flags = procflags.INCONSISTENT_BANDS
1205 flux_flags = [procflags.INCONSISTENT_BANDS] * tot_nband
1206 band_flux = [np.nan] * tot_nband
1207 band_flux_err = [np.nan] * tot_nband
1208 else:
1209 sum_data = _sum_bands_wavg(
1210 all_res=all_res,
1211 all_is_shear_band=all_is_shear_band,
1212 all_wgts=all_wgts,
1213 all_flags=all_flags,
1214 all_wgt_res=None,
1215 )
1216 mdet_flags = copy.copy(sum_data["final_flags"])
1218 psf_sum_data = _sum_bands_wavg(
1219 all_res=all_psf_res,
1220 all_is_shear_band=all_is_shear_band,
1221 all_wgts=all_wgts,
1222 all_flags=all_flags,
1223 all_wgt_res=all_res,
1224 )
1225 psf_flags = copy.copy(psf_sum_data["final_flags"])
1227 if (
1228 sum_data["final_flags"] == 0
1229 and psf_sum_data["final_flags"] == 0
1230 and sum_data["used_shear_bands"] != psf_sum_data["used_shear_bands"]
1231 ):
1232 psf_flags |= procflags.INCONSISTENT_BANDS
1233 mdet_flags |= procflags.INCONSISTENT_BANDS
1235 band_flux = []
1236 band_flux_err = []
1237 flux_flags = []
1238 for gres, flags in zip(all_res, all_flags):
1239 # the input flags mark very basic failures and are ORed across all bands
1240 # these are things like missing and or all zero-weight data, edges, etc.
1241 _flux_flags = 0
1242 _flux_flags |= flags
1244 # an object fit missing in any band is bad too
1245 if gres is None:
1246 _flux_flags |= procflags.MISSING_BAND
1247 band_flux.append(np.nan)
1248 band_flux_err.append(np.nan)
1249 elif (
1250 "flux" not in gres
1251 or "flux_err" not in gres
1252 or "flux_flags" not in gres
1253 ):
1254 _flux_flags |= procflags.NOMOMENTS_FAILURE
1255 band_flux.append(np.nan)
1256 band_flux_err.append(np.nan)
1257 else:
1258 _flux_flags |= gres["flux_flags"]
1259 band_flux.append(gres['flux'])
1260 band_flux_err.append(gres['flux_err'])
1261 flux_flags.append(_flux_flags)
1263 if psf_flags == 0:
1264 psf_raw_mom = psf_sum_data["raw_mom"] / psf_sum_data["wgt_sum"]
1265 psf_raw_mom_cov = psf_sum_data["raw_mom_cov"] / (psf_sum_data["wgt_sum"]**2)
1266 psf_raw_flux = psf_sum_data["flux"] / psf_sum_data["wgt_sum"]
1267 psf_raw_flux_var = psf_sum_data["flux_var"] / (psf_sum_data["wgt_sum"]**2)
1269 psf_momres = _make_mom_res(
1270 raw_mom=psf_raw_mom,
1271 raw_mom_cov=psf_raw_mom_cov,
1272 raw_flux=psf_raw_flux,
1273 raw_flux_var=psf_raw_flux_var,
1274 fwhm_reg=0,
1275 )
1277 psf_flags |= psf_momres["flags"]
1278 data[n("psf_g")] = psf_momres['e']
1279 data[n("psf_T")] = psf_momres['T']
1281 if mdet_flags == 0:
1282 raw_mom = sum_data["raw_mom"] / sum_data["wgt_sum"]
1283 raw_mom_cov = sum_data["raw_mom_cov"] / (sum_data["wgt_sum"]**2)
1284 raw_flux = sum_data["flux"] / sum_data["wgt_sum"]
1285 raw_flux_var = sum_data["flux_var"] / (sum_data["wgt_sum"]**2)
1287 momres = _make_mom_res(
1288 raw_mom=raw_mom,
1289 raw_mom_cov=raw_mom_cov,
1290 raw_flux=raw_flux,
1291 raw_flux_var=raw_flux_var,
1292 fwhm_reg=fwhm_reg,
1293 )
1295 mdet_flags |= momres["flags"]
1296 for col in ['s2n', 'T', 'T_err', 'T_flags']:
1297 data[n(col)] = momres[col]
1298 for col in ['e', 'e_cov']:
1299 data[n(col.replace('e', 'g'))] = momres[col]
1300 if psf_flags == 0:
1301 data[n('T_ratio')] = data[n('T')] / data[n('psf_T')]
1303 if psf_flags != 0:
1304 mdet_flags |= procflags.PSF_FAILURE
1306 if tot_nband > 1:
1307 data[n('band_flux')] = np.array(band_flux)
1308 data[n('band_flux_err')] = np.array(band_flux_err)
1309 data[n('band_flux_flags')] = np.array(flux_flags)
1310 elif tot_nband == 1:
1311 data[n('band_flux')] = band_flux[0]
1312 data[n('band_flux_err')] = band_flux_err[0]
1313 data[n('band_flux_flags')] = flux_flags[0]
1314 else:
1315 data[n('band_flux_flags')] = procflags.MISSING_BAND
1317 # now we set the flags as they would have been set in our moments code
1318 # any PSF failure in a shear band causes a non-zero flags value
1319 data[n('psf_flags')] = psf_flags
1320 data[n('obj_flags')] = mdet_flags
1321 all_flags = mdet_flags
1322 for f in flux_flags:
1323 all_flags |= f
1324 data[n('flags')] = all_flags
1326 if data[n('flags')] != 0:
1327 logger.debug(
1328 "fitter failed: flags = %s",
1329 procflags.get_procflags_str(data[n("flags")])
1330 )
1332 return data
1335def get_wavg_output_struct(nband, model, shear_bands=None):
1336 """
1337 make an output struct with default values set
1339 flags and psf_flags are set to NO_ATTEMPT. The float
1340 fields are set to np.nan
1342 Parameters
1343 ----------
1344 nband: int
1345 Number of bands
1346 model: str
1347 The model or "kind" of fitter
1348 shear_bands : list of int, optional
1349 A list of indices into each mbobs that denotes which band is used for shear.
1350 If given, these are added to the output as a shear_bands field. If not given,
1351 this field is left empty.
1353 Returns
1354 -------
1355 ndarray with fields
1356 """
1357 dt = _make_combine_fit_results_wavg_dtype(
1358 nband=nband, model=model, shear_bands=shear_bands
1359 )
1360 data = np.zeros(1, dtype=dt)
1362 n = Namer(front=model)
1363 data[n('flags')] = procflags.NO_ATTEMPT
1364 data[n('psf_flags')] = procflags.NO_ATTEMPT
1365 data[n("obj_flags")] = procflags.NO_ATTEMPT
1366 data[n("T_flags")] = procflags.NO_ATTEMPT
1367 data[n("band_flux_flags")] = procflags.NO_ATTEMPT
1369 for name in data.dtype.names:
1370 if "flags" not in name:
1371 # all are float except flags
1372 data[name] = np.nan
1374 if shear_bands is not None:
1375 assert len(shear_bands) <= MAX_NUM_SHEAR_BANDS
1376 data["shear_bands"] = "".join("%s" % b for b in sorted(shear_bands))
1378 return data
1381def _make_combine_fit_results_wavg_dtype(nband, model, shear_bands):
1382 n = Namer(front=model)
1383 dt = [
1384 (n("flags"), 'i4'),
1385 (n('psf_flags'), 'i4'),
1386 (n('psf_g'), 'f8', 2),
1387 (n('psf_T'), 'f8'),
1388 (n("obj_flags"), 'i4'),
1389 (n("s2n"), "f8"),
1390 (n("g"), "f8", 2),
1391 (n("g_cov"), "f8", (2, 2)),
1392 (n("T"), "f8"),
1393 (n("T_flags"), "i4"),
1394 (n("T_err"), "f8"),
1395 (n("T_ratio"), "f8"),
1396 ]
1397 if nband > 1:
1398 dt += [
1399 (n("band_flux_flags"), 'i4', nband),
1400 (n("band_flux"), "f8", nband),
1401 (n("band_flux_err"), "f8", nband),
1402 ]
1403 else:
1404 dt += [
1405 (n("band_flux_flags"), 'i4'),
1406 (n("band_flux"), "f8"),
1407 (n("band_flux_err"), "f8"),
1408 ]
1410 dt += [("shear_bands", "U%d" % MAX_NUM_SHEAR_BANDS)]
1412 return dt
1415def symmetrize_obs_weights(obs):
1416 """Applies 4-fold symmetry to zero weight pixels in an observation.
1418 Parameters
1419 ----------
1420 obs : ngmix.Observation
1421 The observation to symmetrize.
1423 Returns
1424 -------
1425 sym_obs : ngmix.Observation
1426 A copy of the input observation with a symmetrized weight map.
1427 """
1428 sym_obs = obs.copy()
1429 if np.any(obs.weight <= 0):
1430 new_wgt = obs.weight.copy()
1431 for k in [1, 2, 3]:
1432 msk = np.rot90(obs.weight, k=k) <= 0
1433 new_wgt[msk] = 0
1435 with sym_obs.writeable():
1436 if not np.any(new_wgt > 0):
1437 sym_obs.ignore_zero_weight = False
1438 sym_obs.weight[:, :] = new_wgt
1440 return sym_obs
1443def combine_fit_res(all_res):
1444 """Combine fit result data structures.
1446 Parameters
1447 ----------
1448 all_res : list of np.ndarray
1449 The list of structured array results.
1451 Returns
1452 -------
1453 res : np.ndarray
1454 The combined results list.
1455 """
1457 if len(all_res) == 1:
1458 return all_res[0]
1460 dupe_cols = [
1461 "shear_bands",
1462 ]
1464 nobj = None
1465 dt = []
1466 for _res in all_res:
1467 if _res is not None:
1468 if nobj is None:
1469 nobj = _res.shape[0]
1470 else:
1471 if nobj != _res.shape[0]:
1472 raise RuntimeError("All fit results must be the same length!")
1474 if len(dt) == 0:
1475 dt.extend(_res.dtype.descr)
1476 else:
1477 for descr in _res.dtype.descr:
1478 if descr[0] not in dupe_cols:
1479 dt.append(descr)
1480 else:
1481 if nobj is None:
1482 nobj = 0
1483 else:
1484 if nobj != 0:
1485 raise RuntimeError(
1486 "All fit results must zero length if one is None!"
1487 )
1489 if nobj > 0:
1490 res = np.zeros(nobj, dtype=dt)
1491 for i, _res in enumerate(all_res):
1492 for col in _res.dtype.names:
1493 if col in dupe_cols:
1494 if i > 0:
1495 if not np.array_equal(res[col], _res[col]):
1496 raise RuntimeError(
1497 "Inconsistent column values "
1498 "for %s when combining results!" % col
1499 )
1500 else:
1501 res[col] = _res[col]
1502 else:
1503 res[col] = _res[col]
1505 return res
1506 else:
1507 if any(_res is not None for _res in all_res):
1508 return np.zeros(0, dtype=dt)
1509 else:
1510 return None