diff --git a/doc/changes.rst b/doc/changes.rst index be1a039d..38a44cec 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -11,6 +11,17 @@ Change Log * Do not constrain the SPS age by default [`PR #132`_]. * Bug fix of emission-line subtracted Dn(4000) measurement [`PR #135`_]. * Update IGM attenuation coefficients [`PR #136`_]. +* Several significant changes [`PR #137`_]: + + * Record the observed-space emission-line amplitude in ``_AMP`` and move the + model-space amplitude to ``_MODELAMP``. + * Demand at least 12 pixels to measure the scatter in the pixels under the + line (therefore ``_AMP_IVAR`` should be more reliable for narrow lines). + * Major bug fix whereby the model emission-line spectra were not being + convolved with the resolution matrix. + * Redefine ``_CHI2`` for an emission line as the observed not reduced chi2. + * Switch from (deprecated) ``pkg_resources`` to ``importlib``. + * Updated documentation (data model) and several non-negligible speed-ups. .. _`PR #115`: https://github.com/desihub/fastspecfit/pull/115 .. _`PR #116`: https://github.com/desihub/fastspecfit/pull/116 @@ -18,6 +29,7 @@ Change Log .. _`PR #132`: https://github.com/desihub/fastspecfit/pull/132 .. _`PR #135`: https://github.com/desihub/fastspecfit/pull/135 .. _`PR #136`: https://github.com/desihub/fastspecfit/pull/136 +.. _`PR #137`: https://github.com/desihub/fastspecfit/pull/137 2.1.2 (2023-04-01) ------------------ diff --git a/doc/fastspec.rst b/doc/fastspec.rst index 024a4206..db6f8fb0 100644 --- a/doc/fastspec.rst +++ b/doc/fastspec.rst @@ -184,6 +184,7 @@ Name Type Units Descript MGII_DOUBLET_RATIO float32 MgII 2796 / 2803 doublet line-ratio. OII_DOUBLET_RATIO float32 [OII] 3726 / 3729 doublet line-ratio. SII_DOUBLET_RATIO float32 [SII] 6731 / 6716 doublet line-ratio. + LYALPHA_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission-line amplitude. LYALPHA_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. LYALPHA_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. LYALPHA_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -198,8 +199,9 @@ Name Type Units Descript LYALPHA_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. LYALPHA_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. LYALPHA_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - LYALPHA_CHI2 float32 Reduced chi^2 of the line-fit. + LYALPHA_CHI2 float32 Chi-squared of the line-fit. LYALPHA_NPIX int32 Number of pixels attributed to the emission line. + OI_1304_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission-line amplitude. OI_1304_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OI_1304_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OI_1304_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -214,8 +216,9 @@ Name Type Units Descript OI_1304_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OI_1304_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OI_1304_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OI_1304_CHI2 float32 Reduced chi^2 of the line-fit. + OI_1304_CHI2 float32 Chi-squared of the line-fit. OI_1304_NPIX int32 Number of pixels attributed to the emission line. + SILIV_1396_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission-line amplitude. SILIV_1396_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SILIV_1396_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SILIV_1396_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -230,8 +233,9 @@ Name Type Units Descript SILIV_1396_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SILIV_1396_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SILIV_1396_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SILIV_1396_CHI2 float32 Reduced chi^2 of the line-fit. + SILIV_1396_CHI2 float32 Chi-squared of the line-fit. SILIV_1396_NPIX int32 Number of pixels attributed to the emission line. + CIV_1549_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. CIV_1549_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. CIV_1549_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. CIV_1549_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -246,8 +250,9 @@ Name Type Units Descript CIV_1549_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. CIV_1549_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. CIV_1549_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - CIV_1549_CHI2 float32 Reduced chi^2 of the line-fit. + CIV_1549_CHI2 float32 Chi-squared of the line-fit. CIV_1549_NPIX int32 Number of pixels attributed to the emission line. + SILIII_1892_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. SILIII_1892_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SILIII_1892_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SILIII_1892_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -262,8 +267,9 @@ Name Type Units Descript SILIII_1892_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SILIII_1892_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SILIII_1892_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SILIII_1892_CHI2 float32 Reduced chi^2 of the line-fit. + SILIII_1892_CHI2 float32 Chi-squared of the line-fit. SILIII_1892_NPIX int32 Number of pixels attributed to the emission line. + CIII_1908_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. CIII_1908_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. CIII_1908_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. CIII_1908_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -278,8 +284,9 @@ Name Type Units Descript CIII_1908_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. CIII_1908_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. CIII_1908_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - CIII_1908_CHI2 float32 Reduced chi^2 of the line-fit. + CIII_1908_CHI2 float32 Chi-squared of the line-fit. CIII_1908_NPIX int32 Number of pixels attributed to the emission line. + MGII_2796_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. MGII_2796_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. MGII_2796_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. MGII_2796_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -293,8 +300,9 @@ Name Type Units Descript MGII_2796_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. MGII_2796_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. MGII_2796_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - MGII_2796_CHI2 float32 Reduced chi^2 of the line-fit. + MGII_2796_CHI2 float32 Chi-squared of the line-fit. MGII_2796_NPIX int32 Number of pixels attributed to the emission line. + MGII_2803_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. MGII_2803_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. MGII_2803_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. MGII_2803_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -309,8 +317,9 @@ Name Type Units Descript MGII_2803_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. MGII_2803_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. MGII_2803_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - MGII_2803_CHI2 float32 Reduced chi^2 of the line-fit. + MGII_2803_CHI2 float32 Chi-squared of the line-fit. MGII_2803_NPIX int32 Number of pixels attributed to the emission line. + NEV_3346_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NEV_3346_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NEV_3346_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NEV_3346_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -325,8 +334,9 @@ Name Type Units Descript NEV_3346_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NEV_3346_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NEV_3346_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NEV_3346_CHI2 float32 Reduced chi^2 of the line-fit. + NEV_3346_CHI2 float32 Chi-squared of the line-fit. NEV_3346_NPIX int32 Number of pixels attributed to the emission line. + NEV_3426_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NEV_3426_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NEV_3426_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NEV_3426_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -341,8 +351,9 @@ Name Type Units Descript NEV_3426_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NEV_3426_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NEV_3426_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NEV_3426_CHI2 float32 Reduced chi^2 of the line-fit. + NEV_3426_CHI2 float32 Chi-squared of the line-fit. NEV_3426_NPIX int32 Number of pixels attributed to the emission line. + OII_3726_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OII_3726_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OII_3726_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OII_3726_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -357,8 +368,9 @@ Name Type Units Descript OII_3726_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OII_3726_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OII_3726_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OII_3726_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + OII_3726_CHI2 float32 Chi-squared of the line-fit. OII_3726_NPIX int32 Number of pixels attributed to the emission line. + OII_3729_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OII_3729_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OII_3729_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OII_3729_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -373,8 +385,9 @@ Name Type Units Descript OII_3729_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OII_3729_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OII_3729_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OII_3729_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + OII_3729_CHI2 float32 Chi-squared of the line-fit. OII_3729_NPIX int32 Number of pixels attributed to the emission line. + NEIII_3869_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NEIII_3869_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NEIII_3869_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NEIII_3869_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -389,8 +402,9 @@ Name Type Units Descript NEIII_3869_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NEIII_3869_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NEIII_3869_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NEIII_3869_CHI2 float32 Reduced chi^2 of the line-fit. + NEIII_3869_CHI2 float32 Chi-squared of the line-fit. NEIII_3869_NPIX int32 Number of pixels attributed to the emission line. + H6_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. H6_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. H6_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. H6_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -405,8 +419,9 @@ Name Type Units Descript H6_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. H6_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. H6_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - H6_CHI2 float32 Reduced chi^2 of the line-fit. + H6_CHI2 float32 Chi-squared of the line-fit. H6_NPIX int32 Number of pixels attributed to the emission line. + H6_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. H6_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. H6_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. H6_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -420,8 +435,9 @@ Name Type Units Descript H6_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. H6_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. H6_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - H6_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + H6_BROAD_CHI2 float32 Chi-squared of the line-fit. H6_BROAD_NPIX int32 Number of pixels attributed to the emission line. + HEPSILON_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HEPSILON_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HEPSILON_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HEPSILON_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -436,8 +452,9 @@ Name Type Units Descript HEPSILON_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HEPSILON_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HEPSILON_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HEPSILON_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HEPSILON_CHI2 float32 Chi-squared of the line-fit. HEPSILON_NPIX int32 Number of pixels attributed to the emission line. + HEPSILON_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HEPSILON_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HEPSILON_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HEPSILON_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -452,8 +469,9 @@ Name Type Units Descript HEPSILON_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HEPSILON_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HEPSILON_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HEPSILON_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HEPSILON_BROAD_CHI2 float32 Chi-squared of the line-fit. HEPSILON_BROAD_NPIX int32 Number of pixels attributed to the emission line. + HDELTA_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HDELTA_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HDELTA_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HDELTA_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -468,8 +486,9 @@ Name Type Units Descript HDELTA_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HDELTA_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HDELTA_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HDELTA_CHI2 float32 Reduced chi^2 of the line-fit. + HDELTA_CHI2 float32 Chi-squared of the line-fit. HDELTA_NPIX int32 Number of pixels attributed to the emission line. + HDELTA_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HDELTA_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HDELTA_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HDELTA_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -484,8 +503,9 @@ Name Type Units Descript HDELTA_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HDELTA_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HDELTA_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HDELTA_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HDELTA_BROAD_CHI2 float32 Chi-squared of the line-fit. HDELTA_BROAD_NPIX int32 Number of pixels attributed to the emission line. + HGAMMA_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HGAMMA_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HGAMMA_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HGAMMA_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -500,8 +520,9 @@ Name Type Units Descript HGAMMA_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HGAMMA_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HGAMMA_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HGAMMA_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HGAMMA_CHI2 float32 Chi-squared of the line-fit. HGAMMA_NPIX int32 Number of pixels attributed to the emission line. + HGAMMA_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HGAMMA_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HGAMMA_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HGAMMA_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -516,8 +537,9 @@ Name Type Units Descript HGAMMA_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HGAMMA_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HGAMMA_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HGAMMA_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HGAMMA_BROAD_CHI2 float32 Chi-squared of the line-fit. HGAMMA_BROAD_NPIX int32 Number of pixels attributed to the emission line. + OIII_4363_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OIII_4363_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OIII_4363_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OIII_4363_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -532,8 +554,9 @@ Name Type Units Descript OIII_4363_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OIII_4363_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OIII_4363_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OIII_4363_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + OIII_4363_CHI2 float32 Chi-squared of the line-fit. OIII_4363_NPIX int32 Number of pixels attributed to the emission line. + HEI_4471_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HEI_4471_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HEI_4471_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HEI_4471_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -548,8 +571,9 @@ Name Type Units Descript HEI_4471_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HEI_4471_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HEI_4471_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HEI_4471_CHI2 float32 Reduced chi^2 of the line-fit. + HEI_4471_CHI2 float32 Chi-squared of the line-fit. HEI_4471_NPIX int32 Number of pixels attributed to the emission line. + HEII_4686_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HEII_4686_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HEII_4686_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HEII_4686_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -564,8 +588,9 @@ Name Type Units Descript HEII_4686_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HEII_4686_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HEII_4686_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HEII_4686_CHI2 float32 Reduced chi^2 of the line-fit. + HEII_4686_CHI2 float32 Chi-squared of the line-fit. HEII_4686_NPIX int32 Number of pixels attributed to the emission line. + HBETA_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HBETA_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HBETA_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HBETA_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -580,8 +605,9 @@ Name Type Units Descript HBETA_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HBETA_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HBETA_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HBETA_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HBETA_CHI2 float32 Chi-squared of the line-fit. HBETA_NPIX int32 Number of pixels attributed to the emission line. + HBETA_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HBETA_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HBETA_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HBETA_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -596,8 +622,9 @@ Name Type Units Descript HBETA_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HBETA_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HBETA_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HBETA_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HBETA_BROAD_CHI2 float32 Chi-squared of the line-fit. HBETA_BROAD_NPIX int32 Number of pixels attributed to the emission line. + OIII_4959_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OIII_4959_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OIII_4959_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OIII_4959_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -612,8 +639,9 @@ Name Type Units Descript OIII_4959_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OIII_4959_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OIII_4959_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OIII_4959_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + OIII_4959_CHI2 float32 Chi-squared of the line-fit. OIII_4959_NPIX int32 Number of pixels attributed to the emission line. + OIII_5007_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OIII_5007_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OIII_5007_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OIII_5007_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -628,8 +656,9 @@ Name Type Units Descript OIII_5007_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OIII_5007_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OIII_5007_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OIII_5007_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + OIII_5007_CHI2 float32 Chi-squared of the line-fit. OIII_5007_NPIX int32 Number of pixels attributed to the emission line. + NII_5755_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NII_5755_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NII_5755_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NII_5755_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -644,8 +673,9 @@ Name Type Units Descript NII_5755_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NII_5755_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NII_5755_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NII_5755_CHI2 float32 Reduced chi^2 of the line-fit. + NII_5755_CHI2 float32 Chi-squared of the line-fit. NII_5755_NPIX int32 Number of pixels attributed to the emission line. + HEI_5876_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HEI_5876_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HEI_5876_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HEI_5876_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -660,8 +690,9 @@ Name Type Units Descript HEI_5876_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HEI_5876_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HEI_5876_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HEI_5876_CHI2 float32 Reduced chi^2 of the line-fit. + HEI_5876_CHI2 float32 Chi-squared of the line-fit. HEI_5876_NPIX int32 Number of pixels attributed to the emission line. + OI_6300_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OI_6300_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OI_6300_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OI_6300_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -676,8 +707,9 @@ Name Type Units Descript OI_6300_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OI_6300_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OI_6300_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OI_6300_CHI2 float32 Reduced chi^2 of the line-fit. + OI_6300_CHI2 float32 Chi-squared of the line-fit. OI_6300_NPIX int32 Number of pixels attributed to the emission line. + NII_6548_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NII_6548_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NII_6548_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NII_6548_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -692,8 +724,9 @@ Name Type Units Descript NII_6548_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NII_6548_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NII_6548_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NII_6548_CHI2 float32 Reduced chi^2 of the line-fit. + NII_6548_CHI2 float32 Chi-squared of the line-fit. NII_6548_NPIX int32 Number of pixels attributed to the emission line. + HALPHA_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HALPHA_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HALPHA_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HALPHA_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -708,8 +741,9 @@ Name Type Units Descript HALPHA_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HALPHA_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HALPHA_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HALPHA_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HALPHA_CHI2 float32 Chi-squared of the line-fit. HALPHA_NPIX int32 Number of pixels attributed to the emission line. + HALPHA_BROAD_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. HALPHA_BROAD_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. HALPHA_BROAD_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. HALPHA_BROAD_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -724,8 +758,9 @@ Name Type Units Descript HALPHA_BROAD_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. HALPHA_BROAD_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. HALPHA_BROAD_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - HALPHA_BROAD_CHI2 float32 Reduced chi^2 of the line-fit (default value 1e6). + HALPHA_BROAD_CHI2 float32 Chi-squared of the line-fit. HALPHA_BROAD_NPIX int32 Number of pixels attributed to the emission line. + NII_6584_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. NII_6584_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. NII_6584_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. NII_6584_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -740,8 +775,9 @@ Name Type Units Descript NII_6584_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. NII_6584_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. NII_6584_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - NII_6584_CHI2 float32 Reduced chi^2 of the line-fit. + NII_6584_CHI2 float32 Chi-squared of the line-fit. NII_6584_NPIX int32 Number of pixels attributed to the emission line. + SII_6716_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. SII_6716_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SII_6716_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SII_6716_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -756,8 +792,9 @@ Name Type Units Descript SII_6716_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SII_6716_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SII_6716_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SII_6716_CHI2 float32 Reduced chi^2 of the line-fit. + SII_6716_CHI2 float32 Chi-squared of the line-fit. SII_6716_NPIX int32 Number of pixels attributed to the emission line. + SII_6731_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. SII_6731_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SII_6731_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SII_6731_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -772,8 +809,9 @@ Name Type Units Descript SII_6731_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SII_6731_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SII_6731_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SII_6731_CHI2 float32 Reduced chi^2 of the line-fit. + SII_6731_CHI2 float32 Chi-squared of the line-fit. SII_6731_NPIX int32 Number of pixels attributed to the emission line. + OII_7320_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OII_7320_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OII_7320_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OII_7320_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -788,8 +826,9 @@ Name Type Units Descript OII_7320_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OII_7320_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OII_7320_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OII_7320_CHI2 float32 Reduced chi^2 of the line-fit. + OII_7320_CHI2 float32 Chi-squared of the line-fit. OII_7320_NPIX int32 Number of pixels attributed to the emission line. + OII_7330_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. OII_7330_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. OII_7330_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. OII_7330_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -804,8 +843,9 @@ Name Type Units Descript OII_7330_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. OII_7330_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. OII_7330_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - OII_7330_CHI2 float32 Reduced chi^2 of the line-fit. + OII_7330_CHI2 float32 Chi-squared of the line-fit. OII_7330_NPIX int32 Number of pixels attributed to the emission line. + SIII_9069_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. SIII_9069_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SIII_9069_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SIII_9069_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -820,8 +860,9 @@ Name Type Units Descript SIII_9069_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SIII_9069_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SIII_9069_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SIII_9069_CHI2 float32 Reduced chi^2 of the line-fit. + SIII_9069_CHI2 float32 Chi-squared of the line-fit. SIII_9069_NPIX int32 Number of pixels attributed to the emission line. + SIII_9532_MODELAMP float32 1e-17 erg / (Angstrom cm2 s) Model emission line amplitude. SIII_9532_AMP float32 1e-17 erg / (Angstrom cm2 s) Emission line amplitude. SIII_9532_AMP_IVAR float32 1e+34 Angstrom2 cm4 s2 / erg2 Inverse variance of line-amplitude. SIII_9532_FLUX float32 1e-17 erg / (cm2 s) Gaussian-integrated emission-line flux. @@ -836,7 +877,7 @@ Name Type Units Descript SIII_9532_EW_IVAR float32 1 / Angstrom2 Inverse variance of equivalent width. SIII_9532_FLUX_LIMIT float32 erg / (cm2 s) One-sigma upper limit on the emission line flux. SIII_9532_EW_LIMIT float32 Angstrom One-sigma upper limit on the emission line equivalent width. - SIII_9532_CHI2 float32 Reduced chi^2 of the line-fit. + SIII_9532_CHI2 float32 Chi-squared of the line-fit. SIII_9532_NPIX int32 Number of pixels attributed to the emission line. ============================ ============ ============================= ============================================ diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 63efe732..f4e6b505 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -14,61 +14,6 @@ from fastspecfit.io import FLUXNORM from fastspecfit.util import C_LIGHT -#import numba -#@numba.jit(nopython=True) -def nnls(A, b, maxiter=None, eps=1e-7, v=False): - # https://en.wikipedia.org/wiki/Non-negative_least_squares#Algorithms - # not sure what eps should be set to - m, n = A.shape - if b.shape[0] != m: - raise ValueError()#f"Shape mismatch: {A.shape} {b.shape}. Expected: (m, n) (m) ") - if maxiter is None: - # this is what scipy does when maxiter is not specified - maxiter = 3 * n - # init - P = np.zeros(n).astype(bool) - R = np.ones(n).astype(bool) - x = np.zeros(n) - s = np.zeros(n) - # compute these once ahead of time - ATA = A.T.dot(A) - ATb = A.T.dot(b) - # main loop - w = ATb - ATA.dot(x) - j = np.argmax(R * w) - c = 0 # iteration count - # while R != {} and max(w[R]) > eps - while np.any(R) and w[j] > eps: - #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}") - # add j to P, remove j from R - P[j], R[j] = True, False - s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P]) - s[R] = 0 - d = 0 # inner loop iteration count, for debugging - #if v: print(f"{c=}", f"{P=}\n {R=}\n {s=}") - # while P != {} and min(s[P]) < eps - # make sure P is not empty before checking min s[P] - while np.any(P) and np.min(s[P]) < eps: - i = P & (s < eps) - #if v: print(f" {d=}", f" {P=}\n {i=}\n {s=}\n {x=}") - a = np.min(x[i] / (x[i] - s[i])) - x = x + a * (s - x) - j = P & (x < eps) - R[j], P[j] = True, False - s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P]) - s[R] = 0 - d += 1 - x[:] = s - # w = A.T.dot(b - A.dot(x)) - w = ATb - ATA.dot(x) - j = np.argmax(R * w) - #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}") - c += 1 - if c >= maxiter: - break - res = np.linalg.norm(A.dot(x) - b) - return x, res - def _smooth_continuum(wave, flux, ivar, redshift, medbin=150, smooth_window=50, smooth_step=10, maskkms_uv=3000.0, maskkms_balmer=1000.0, maskkms_narrow=200.0, @@ -1485,7 +1430,7 @@ def templates2data(self, _templateflux, _templatewave, redshift=0.0, dluminosity @staticmethod def call_nnls(modelflux, flux, ivar, xparam=None, debug=False, - interpolate_coeff=False, xlabel=None, log=None): + interpolate_coeff=False, xlabel=None, png=None, log=None): """Wrapper on nnls. Works with both spectroscopic and photometric input and with both 2D and @@ -1535,11 +1480,14 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False, try: imin = find_minima(chi2grid)[0] - xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2]) + if debug: + xbest, xerr, chi2min, warn, (a, b, c) = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], return_coeff=True) + else: + xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2]) except: errmsg = 'A problem was encountered minimizing chi2.' log.warning(errmsg) - imin, xbest, xerr, chi2min, warn = 0, 0.0, 0.0, 0.0, 1 + imin, xbest, xerr, chi2min, warn, (a, b, c) = 0, 0.0, 0.0, 0.0, 1, (0., 0., 0.) if warn == 0: xivar = 1.0 / xerr**2 @@ -1565,29 +1513,37 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False, if debug: if xivar > 0: - leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min) + leg = r'${:.1f}\pm{:.1f}$'.format(xbest, 1 / np.sqrt(xivar)) + #leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min) else: leg = r'${:.3f}$'.format(xbest) + if xlabel: + leg = f'{xlabel} = {leg}' + import matplotlib.pyplot as plt + import seaborn as sns + sns.set(context='talk', style='ticks', font_scale=0.8) fig, ax = plt.subplots() - ax.scatter(xparam, chi2grid) - ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], color='red') - #ax.set_ylim(chi2min*0.95, np.max(chi2grid[imin-1:imin+2])*1.05) - #ax.plot(xx, np.polyval([aa, bb, cc], xx), ls='--') - ax.axvline(x=xbest, color='k') - if xivar > 0: - ax.axhline(y=chi2min, color='k') - #ax.set_yscale('log') - #ax.set_ylim(chi2min, 63.3) + ax.scatter(xparam, chi2grid-chi2min, marker='s', s=50, color='gray', edgecolor='k') + #ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2]-chi2min, color='red') + if not np.all(np.array([a, b, c]) == 0): + ax.plot(xparam, np.polyval([a, b, c], xparam)-chi2min, lw=2, ls='--') + #ax.axvline(x=xbest, color='k') + #if xivar > 0: + # ax.axhline(y=chi2min, color='k') if xlabel: ax.set_xlabel(xlabel) #ax.text(0.03, 0.9, '{}={}'.format(xlabel, leg), ha='left', # va='center', transform=ax.transAxes) - ax.text(0.03, 0.9, leg, ha='left', va='center', transform=ax.transAxes) - ax.set_ylabel(r'$\chi^2$') + ax.text(0.9, 0.9, leg, ha='right', va='center', transform=ax.transAxes) + ax.set_ylabel(r'$\Delta\chi^2$') + #ax.set_ylabel(r'$\chi^2 - {:.1f}$'.format(chi2min)) #fig.savefig('qa-chi2min.png') - fig.savefig('desi-users/ioannis/tmp/qa-chi2min.png') + fig.tight_layout() + if png: + log.info(f'Writing {png}') + fig.savefig(png) return chi2min, xbest, xivar, bestcoeff @@ -1797,7 +1753,7 @@ def _kcorr_and_absmag(filters_out, band_shift): return kcorr, absmag, ivarabsmag, bestmaggies, lums, cfluxes def continuum_specfit(data, result, templatecache, nophoto=False, constrain_age=False, - fastphot=False, log=None, verbose=False): + fastphot=False, log=None, verbose=False, debug_plots=False): """Fit the non-negative stellar continuum of a single spectrum. Parameters @@ -1960,9 +1916,10 @@ def _younger_than_universe(age, tuniv, agepad=0.5): vdispchi2min, vdispbest, vdispivar, _ = CTools.call_nnls( ztemplateflux_vdisp[Ivdisp, :, :], specflux[Ivdisp], specivar[Ivdisp], - xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', debug=False, log=log) + xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', log=log, + debug=debug_plots, png='deltachi2-vdisp.png') log.info('Fitting for the velocity dispersion took {:.2f} seconds.'.format(time.time()-t0)) - + if vdispivar > 0: # Require vdisp to be measured with S/N>1, which protects # against tiny ivar becomming infinite in the output table. diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 4fe38619..56104eef 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -19,16 +19,15 @@ def read_emlines(): """Read the set of emission lines of interest. """ - from pkg_resources import resource_filename - - linefile = resource_filename('fastspecfit', 'data/emlines.ecsv') + from importlib import resources + linefile = resources.files('fastspecfit').joinpath('data/emlines.ecsv') linetable = Table.read(linefile, format='ascii.ecsv', guess=False) return linetable #import numba #@numba.jit(nopython=True) -def build_emline_model(log10wave, redshift, lineamps, linevshifts, linesigmas, +def build_emline_model(dlog10wave, redshift, lineamps, linevshifts, linesigmas, linewaves, emlinewave, resolution_matrix, camerapix=None): """Given parameters, build the model emission-line spectrum. @@ -37,7 +36,7 @@ def build_emline_model(log10wave, redshift, lineamps, linevshifts, linesigmas, """ from fastspecfit.util import trapz_rebin - log10model = np.zeros_like(log10wave) # [erg/s/cm2/A, observed frame] + #log10model = np.zeros_like(log10wave) # [erg/s/cm2/A, observed frame] # Cut to lines with non-zero amplitudes. #I = linesigmas > 0 @@ -48,15 +47,38 @@ def build_emline_model(log10wave, redshift, lineamps, linevshifts, linesigmas, lineamps = lineamps[I] linewaves = linewaves[I] + # demand at least 20 km/s for rendering the model + if np.any(linesigmas) < 20.: + linesigmas[linesigmas<20.] = 20. + # line-width [log-10 Angstrom] and redshifted wavelength [log-10 Angstrom] log10sigmas = linesigmas / C_LIGHT / np.log(10) linezwaves = np.log10(linewaves * (1.0 + redshift + linevshifts / C_LIGHT)) + # Build the wavelength vector on-the-fly. + if camerapix is None: + minwave = emlinewave[0][0]-2. + maxwave = emlinewave[-1][-1]+2. + else: + minwave = emlinewave[0]-2. + maxwave = emlinewave[-1]+2. + + log10wave = [] + for linezwave, log10sigma in zip(linezwaves, log10sigmas): + log10wave.append(np.arange(linezwave - (5 * log10sigma), linezwave + (5 * log10sigma), dlog10wave)) + log10wave = np.hstack([np.log10(minwave), np.log10(maxwave), ] + log10wave) + S = np.argsort(log10wave) + log10wave = log10wave[S] + log10model = np.zeros_like(log10wave) + for lineamp, linezwave, log10sigma in zip(lineamps, linezwaves, log10sigmas): - J = np.abs(log10wave - linezwave) < (8 * log10sigma) # cut to pixels within +/-N-sigma - if np.count_nonzero(J) > 0: - #print(lineamp, 10**linezwave, 10**log10wave[J].min(), 10**log10wave[J].max()) - log10model[J] = log10model[J] + lineamp * np.exp(-0.5 * (log10wave[J]-linezwave)**2 / log10sigma**2) + J = np.abs(log10wave - linezwave) < (5 * log10sigma) + #print(lineamp, 10**linezwave, 10**log10wave[J].min(), 10**log10wave[J].max()) + log10model[J] += lineamp * np.exp(-0.5 * (log10wave[J]-linezwave)**2 / log10sigma**2) + + #import matplotlib.pyplot as plt + #plt.clf() ; plt.scatter(10**log10wave, log10model, s=1) ; plt.xlim(7500, 8100) ; plt.savefig('desi-users/ioannis/tmp/junk.png') + #pdb.set_trace() # Optionally split into cameras, resample, and convolve with the resolution # matrix. @@ -64,18 +86,18 @@ def build_emline_model(log10wave, redshift, lineamps, linevshifts, linesigmas, if camerapix is None: for icam, specwave in enumerate(emlinewave): _emlinemodel = trapz_rebin(10**log10wave, log10model, specwave) - _emlinemomdel = resolution_matrix[icam].dot(_emlinemodel) + _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) emlinemodel.append(_emlinemodel) return emlinemodel else: for icam, campix in enumerate(camerapix): _emlinemodel = trapz_rebin(10**log10wave, log10model, emlinewave[campix[0]:campix[1]]) - _emlinemomdel = resolution_matrix[icam].dot(_emlinemodel) + _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) emlinemodel.append(_emlinemodel) return np.hstack(emlinemodel) def _objective_function(free_parameters, emlinewave, emlineflux, weights, redshift, - log10wave, resolution_matrix, camerapix, parameters, Ifree, + dlog10wave, resolution_matrix, camerapix, parameters, Ifree, Itied, tiedtoparam, tiedfactor, doubletindx, doubletpair, linewaves): """The parameters array should only contain free (not tied or fixed) parameters.""" @@ -104,7 +126,7 @@ def _objective_function(free_parameters, emlinewave, emlineflux, weights, redshi lineamps[doubletindx] *= lineamps[doubletpair] # Build the emission-line model. - emlinemodel = build_emline_model(log10wave, redshift, lineamps, linevshifts, + emlinemodel = build_emline_model(dlog10wave, redshift, lineamps, linevshifts, linesigmas, linewaves, emlinewave, resolution_matrix, camerapix) @@ -158,9 +180,8 @@ def __init__(self, minspecwave=3500.0, maxspecwave=9900.0, targetid=None): self.linetable = read_emlines() - self.emwave_pixkms = 5.0 # pixel size for internal wavelength array [km/s] - self.dlogwave = self.emwave_pixkms / C_LIGHT / np.log(10) # pixel size [log-lambda] - self.log10wave = np.arange(np.log10(minspecwave), np.log10(maxspecwave), self.dlogwave) + self.emwave_pixkms = 5. # pixel size for internal wavelength array [km/s] + self.dlog10wave = self.emwave_pixkms / C_LIGHT / np.log(10) # pixel size [log-lambda] # default line-sigma for computing upper limits self.limitsigma_narrow = 75.0 @@ -352,6 +373,7 @@ def _fix_parameters(linemodel, verbose=False): final_linemodel['bounds'] = np.zeros((nparam, 2), 'f8') final_linemodel['initial'] = np.zeros(nparam, 'f8') final_linemodel['value'] = np.zeros(nparam, 'f8') + final_linemodel['obsvalue'] = np.zeros(nparam, 'f8') final_linemodel['civar'] = np.zeros(nparam, 'f8') # continuum inverse variance final_linemodel['doubletpair'][self.doubletindx] = self.doubletpair @@ -788,12 +810,13 @@ def populate_linemodel(linemodel, initial_guesses, param_bounds, log, broadbalme linemodel['value'] = linemodel['initial'] # copy def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, - resolution_matrix, camerapix, log=None, verbose=False, debug=False): + resolution_matrix, camerapix, log=None, get_finalamp=False, + verbose=False, debug=False): """Wrapper to call the least-squares minimization given a linemodel. """ from scipy.optimize import least_squares - from fastspecfit.emlines import _objective_function + from fastspecfit.util import trapz_rebin if log is None: from desiutil.log import get_logger, DEBUG @@ -806,7 +829,7 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, linewaves) = self._linemodel_to_parameters(linemodel, self.fit_linetable) log.debug('Optimizing {} free parameters'.format(len(Ifree))) - farg = (emlinewave, emlineflux, weights, redshift, self.log10wave, + farg = (emlinewave, emlineflux, weights, redshift, self.dlog10wave, resolution_matrix, camerapix, parameters, ) + \ (Ifree, Itied, tiedtoparam, tiedfactor, doubletindx, doubletpair, linewaves) @@ -847,9 +870,10 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, # line-amplitude is dropped, too (see MgII 2796 on # sv1-bright-17680-39627622543528153). drop2 = np.zeros(len(parameters), bool) - drop2[Ifree] = parameters[Ifree] == linemodel['value'][Ifree] - drop2 *= notfixed - + if False: + drop2[Ifree] = parameters[Ifree] == linemodel['value'][Ifree] + drop2 *= notfixed + sigmadropped = np.where(self.sigma_param_bool * drop2)[0] if len(sigmadropped) > 0: for lineindx, dropline in zip(sigmadropped, linemodel[sigmadropped]['linename']): @@ -872,9 +896,6 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, log.debug('Dropping {} parameters which are out-of-bounds.'.format(np.sum(drop3))) Idrop = np.where(np.logical_or.reduce((drop1, drop2, drop3)))[0] - #if debug: - # pdb.set_trace() - ## If we are fitting the broad Balmer lines, do some additional sanity checks. if len(Idrop) > 0: log.debug(' Dropping {} unique parameters.'.format(len(Idrop))) @@ -896,17 +917,63 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, out_linemodel['value'] = parameters out_linemodel.meta['nfev'] = fit_info['nfev'] - if False: - bestfit = self.bestfit(out_linemodel, redshift, emlinewave, resolution_matrix, camerapix) - import matplotlib.pyplot as plt - plt.clf() - plt.plot(emlinewave, emlineflux) - plt.plot(emlinewave, bestfit) - #plt.xlim(5800, 6200) - #plt.xlim(6600, 6950) - plt.xlim(5050, 5120) - #plt.xlim(8850, 9050) - plt.savefig('junk.png') + # Get the final line-amplitudes, after resampling and convolution (see + # https://github.com/desihub/fastspecfit/issues/139). Some repeated code + # from build_emline_model... + if get_finalamp: + lineamps, linevshifts, linesigmas = np.array_split(parameters, 3) # 3 parameters per line + lineindxs = np.arange(len(lineamps)) + + I = lineamps > 0 + if np.count_nonzero(I) > 0: + linevshifts = linevshifts[I] + linesigmas = linesigmas[I] + lineamps = lineamps[I] + linewaves = linewaves[I] + lineindxs = lineindxs[I] + + # demand at least 20 km/s for rendering the model + if np.any(linesigmas) < 20.: + linesigmas[linesigmas<20.] = 20. + + if camerapix is None: + minwave = emlinewave[0][0]-2. + maxwave = emlinewave[-1][-1]+2. + else: + minwave = emlinewave[0]-2. + maxwave = emlinewave[-1]+2. + + _emlinewave = [] + for icam, campix in enumerate(camerapix): + _emlinewave.append(emlinewave[campix[0]:campix[1]]) + + # line-width [log-10 Angstrom] and redshifted wavelength [log-10 Angstrom] + log10sigmas = linesigmas / C_LIGHT / np.log(10) + linezwaves = np.log10(linewaves * (1.0 + redshift + linevshifts / C_LIGHT)) + + for lineindx, lineamp, linezwave, log10sigma in zip(lineindxs, lineamps, linezwaves, log10sigmas): + log10wave = np.arange(linezwave - (5 * log10sigma), linezwave + (5 * log10sigma), self.dlog10wave) + log10wave = np.hstack((np.log10(minwave), log10wave, np.log10(maxwave))) + log10model = lineamp * np.exp(-0.5 * (log10wave-linezwave)**2 / log10sigma**2) + # Determine which camera we're on and then resample and + # convolve with the resolution matrix. + icam = np.argmin([np.abs((np.max(emwave)-np.min(emwave))/2+np.min(emwave)-10**linezwave) for emwave in _emlinewave]) + model_resamp = trapz_rebin(10**log10wave, log10model, _emlinewave[icam]) + model_convol = resolution_matrix[icam].dot(model_resamp) + out_linemodel['obsvalue'][lineindx] = np.max(model_convol) + #if out_linemodel[lineindx]['param_name'] == 'oiii_5007_amp': + # import matplotlib.pyplot as plt + # bestfit = self.bestfit(out_linemodel, redshift, emlinewave, resolution_matrix, camerapix) + # plt.clf() + # plt.plot(emlinewave, emlineflux, label='data', color='gray', lw=4) + # plt.plot(emlinewave, bestfit, label='bestfit', ls='--', lw=3, alpha=0.7, color='k') + # plt.plot(10**self.log10wave, log10model, label='hires model') + # plt.plot(_emlinewave[icam], model_resamp, label='resamp') + # plt.plot(_emlinewave[icam], model_convol, label='convol', lw=2) + # plt.xlim(5386, 5394) + # plt.legend() + # plt.savefig('desi-users/ioannis/tmp/junk.png') + # pdb.set_trace() return out_linemodel @@ -945,10 +1012,8 @@ def bestfit(self, linemodel, redshift, emlinewave, resolution_matrix, camerapix) # doublets lineamps[doubletindx] *= lineamps[doubletpair] - emlinemodel = build_emline_model(self.log10wave, redshift, lineamps, - linevshifts, linesigmas, linewaves, - emlinewave, resolution_matrix, - camerapix) + emlinemodel = build_emline_model(self.dlog10wave, redshift, lineamps, linevshifts, linesigmas, + linewaves, emlinewave, resolution_matrix, camerapix) return emlinemodel @@ -964,7 +1029,12 @@ def emlinemodel_bestfit(self, specwave, specres, fastspecfit_table, redshift=Non linewaves = self.linetable['restwave'].data - parameters = [fastspecfit_table[param.upper()] for param in self.param_names] + parameters = [] + for param in self.param_names: + if '_amp' in param: + param = param.replace('_amp', '_modelamp') + parameters.append(fastspecfit_table[param.upper()]) + #parameters = [fastspecfit_table[param.upper()] for param in self.param_names] lineamps, linevshifts, linesigmas = np.array_split(parameters, 3) # 3 parameters per line @@ -972,7 +1042,7 @@ def emlinemodel_bestfit(self, specwave, specres, fastspecfit_table, redshift=Non # amplitude parameters are always in the first third of parameters. lineamps[self.doubletindx] *= lineamps[self.doubletpair] - emlinemodel = build_emline_model(self.log10wave, redshift, lineamps, + emlinemodel = build_emline_model(self.dlog10wave, redshift, lineamps, linevshifts, linesigmas, linewaves, specwave, specres, None) @@ -980,7 +1050,7 @@ def emlinemodel_bestfit(self, specwave, specres, fastspecfit_table, redshift=Non def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, emlineivar, oemlineivar, specflux_nolines, redshift, - resolution_matrix, camerapix, log): + resolution_matrix, camerapix, log, nminpix=7): """Populate the output table with the emission-line measurements. """ @@ -989,19 +1059,29 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, for param in finalfit: val = param['value'] + obsval = param['obsvalue'] # special case the tied doublets if param['param_name'] == 'oii_doublet_ratio': result['OII_DOUBLET_RATIO'] = val - result['OII_3726_AMP'] = val * finalfit[param['doubletpair']]['value'] + result['OII_3726_MODELAMP'] = val * finalfit[param['doubletpair']]['value'] + result['OII_3726_AMP'] = val * finalfit[param['doubletpair']]['obsvalue'] elif param['param_name'] == 'sii_doublet_ratio': result['SII_DOUBLET_RATIO'] = val - result['SII_6731_AMP'] = val * finalfit[param['doubletpair']]['value'] + result['SII_6731_MODELAMP'] = val * finalfit[param['doubletpair']]['value'] + result['SII_6731_AMP'] = val * finalfit[param['doubletpair']]['obsvalue'] elif param['param_name'] == 'mgii_doublet_ratio': result['MGII_DOUBLET_RATIO'] = val - result['MGII_2796_AMP'] = val * finalfit[param['doubletpair']]['value'] + result['MGII_2796_MODELAMP'] = val * finalfit[param['doubletpair']]['value'] + result['MGII_2796_AMP'] = val * finalfit[param['doubletpair']]['obsvalue'] else: - result[param['param_name'].upper()] = val + if '_amp' in param['param_name']: + result[param['param_name'].upper().replace('AMP', 'MODELAMP')] = val + result[param['param_name'].upper()] = obsval + else: + result[param['param_name'].upper()] = val + dpixwave = emlinewave[1]-emlinewave[0] # pixel size [Angstrom] + # get continuum fluxes, EWs, and upper limits narrow_sigmas, broad_sigmas, uv_sigmas = [], [], [] narrow_redshifts, broad_redshifts, uv_redshifts = [], [], [] @@ -1023,55 +1103,42 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, linesigma = self.limitsigma_broad linesigma_ang = linesigma * linezwave / C_LIGHT # [observed-frame Angstrom] - #log10sigma = linesigma / C_LIGHT / np.log(10) # line-width [log-10 Angstrom] + + # require at least 3 pixels + if linesigma_ang < 2 * dpixwave: + linesigma_ang_window = 2 * dpixwave + else: + linesigma_ang_window = linesigma_ang # Are the pixels based on the original inverse spectrum fully # masked? If so, set everything to zero and move onto the next line. - lineindx = np.where((emlinewave >= (linezwave - 3.0*linesigma_ang)) * - (emlinewave <= (linezwave + 3.0*linesigma_ang)))[0] + lineindx = np.where((emlinewave >= (linezwave - 3.0*linesigma_ang_window)) * + (emlinewave <= (linezwave + 3.0*linesigma_ang_window)))[0] if len(lineindx) > 0 and np.sum(oemlineivar[lineindx] == 0) / len(lineindx) > 0.3: # use original ivar result['{}_AMP'.format(linename)] = 0.0 + result['{}_MODELAMP'.format(linename)] = 0.0 result['{}_VSHIFT'.format(linename)] = 0.0 result['{}_SIGMA'.format(linename)] = 0.0 else: # number of pixels, chi2, and boxcar integration - lineindx = np.where((emlinewave >= (linezwave - 3.0*linesigma_ang)) * - (emlinewave <= (linezwave + 3.0*linesigma_ang)) * + lineindx = np.where((emlinewave >= (linezwave - 3.0*linesigma_ang_window)) * + (emlinewave <= (linezwave + 3.0*linesigma_ang_window)) * (emlineivar > 0))[0] - - # can happen if sigma is very small (depending on the wavelength) - #if (linezwave > np.min(emlinewave)) * (linezwave < np.max(emlinewave)) * len(lineindx) > 0 and len(lineindx) <= 3: - if (linezwave > np.min(emlinewave)) * (linezwave < np.max(emlinewave)) * (len(lineindx) <= 3): - dwave = emlinewave - linezwave - lineindx = np.argmin(np.abs(dwave)) - if dwave[lineindx] > 0: - pad = np.array([-2, -1, 0, +1]) - else: - pad = np.array([-1, 0, +1, +2]) - - # check to make sure we don't hit the edges - if (lineindx-pad[0]) < 0 or (lineindx+pad[-1]) >= len(emlineivar): - lineindx = np.array([]) - else: - lineindx += pad - # the padded pixels can have ivar==0 - good = oemlineivar[lineindx] > 0 # use the original ivar - lineindx = lineindx[good] - + npix = len(lineindx) result['{}_NPIX'.format(linename)] = npix - if npix > 3: # magic number: required at least XX unmasked pixels centered on the line - + if npix >= nminpix: # magic number: required at least XX unmasked pixels centered on the line + if np.any(emlineivar[lineindx] == 0): errmsg = 'Ivar should never be zero within an emission line!' log.critical(errmsg) raise ValueError(errmsg) # boxcar integration of the flux; should we weight by the line-profile??? - boxflux = np.sum(emlineflux[lineindx]) - boxflux_ivar = 1 / np.sum(1 / emlineivar[lineindx]) + boxflux = np.trapz(emlineflux[lineindx], x=emlinewave[lineindx]) + boxflux_ivar = 1 / np.trapz(1 / emlineivar[lineindx], x=emlinewave[lineindx]) result['{}_BOXFLUX'.format(linename)] = boxflux # * u.erg/(u.second*u.cm**2) result['{}_BOXFLUX_IVAR'.format(linename)] = boxflux_ivar # * u.second**2*u.cm**4/u.erg**2 @@ -1085,48 +1152,30 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, if amp_sigma > 0: result['{}_AMP_IVAR'.format(linename)] = 1 / amp_sigma**2 # * u.second**2*u.cm**4*u.Angstrom**2/u.erg**2 - #if np.isinf(result['{}_AMP_IVAR'.format(linename)]): - # pdb.set_trace() - # require amp > 0 (line not dropped) to compute the flux and chi2 - if result['{}_AMP'.format(linename)] > 0: + if result['{}_MODELAMP'.format(linename)] > 0: # get the emission-line flux linenorm = np.sqrt(2.0 * np.pi) * linesigma_ang # * u.Angstrom - result['{}_FLUX'.format(linename)] = result['{}_AMP'.format(linename)] * linenorm - - #result['{}_FLUX_IVAR'.format(linename)] = result['{}_AMP_IVAR'.format(linename)] / linenorm**2 - #weight = np.exp(-0.5 * np.log10(emlinewave/linezwave)**2 / log10sigma**2) - #weight = (weight / np.max(weight)) > 1e-3 - #result['{}_FLUX_IVAR'.format(linename)] = 1 / np.sum(1 / emlineivar[weight]) - #result['{}_FLUX_IVAR'.format(linename)] = boxflux_ivar # * u.second**2*u.cm**4/u.erg**2 - - # weight by the per-pixel inverse variance line-profile - lineprofile = build_emline_model(self.log10wave, redshift, np.array([result['{}_AMP'.format(linename)]]), - np.array([result['{}_VSHIFT'.format(linename)]]), np.array([result['{}_SIGMA'.format(linename)]]), - np.array([oneline['restwave']]), emlinewave, resolution_matrix, camerapix) + result['{}_FLUX'.format(linename)] = result['{}_MODELAMP'.format(linename)] * linenorm + + ## weight by the per-pixel inverse variance line-profile + #lineprofile = build_emline_model(self.dlog10wave, redshift, np.array([result['{}_MODELAMP'.format(linename)]]), + # np.array([result['{}_VSHIFT'.format(linename)]]), np.array([result['{}_SIGMA'.format(linename)]]), + # np.array([oneline['restwave']]), emlinewave, resolution_matrix, camerapix) + # + #weight = np.sum(lineprofile[lineindx]) + #if weight == 0.0: + # errmsg = 'Line-profile should never sum to zero!' + # log.critical(errmsg) + # raise ValueError(errmsg) + # + #flux_ivar = weight / np.sum(lineprofile[lineindx] / emlineivar[lineindx]) + #result['{}_FLUX_IVAR'.format(linename)] = flux_ivar # * u.second**2*u.cm**4/u.erg**2 - weight = np.sum(lineprofile[lineindx]) - if weight == 0.0: - errmsg = 'Line-profile should never sum to zero!' - log.critical(errmsg) - raise ValueError(errmsg) - - flux_ivar = weight / np.sum(lineprofile[lineindx] / emlineivar[lineindx]) - result['{}_FLUX_IVAR'.format(linename)] = flux_ivar # * u.second**2*u.cm**4/u.erg**2 - - ##if linename == 'OII_3729': - #_indx = np.arange((lineindx[-1]+20)-(lineindx[0]-20))+(lineindx[0]-20) - #import matplotlib.pyplot as plt - #plt.clf() - #plt.plot(emlinewave[_indx], emlineflux[_indx], color='gray') - #plt.plot(emlinewave[_indx], lineprofile[_indx], color='red') - #plt.savefig('desi-users/ioannis/tmp/junk-{}.png'.format(linename)) - - dof = npix - 3 # ??? [redshift, sigma, and amplitude] - chi2 = np.sum(emlineivar[lineindx]*(emlineflux[lineindx]-finalmodel[lineindx])**2) / dof - - result['{}_CHI2'.format(linename)] = chi2 + result['{}_FLUX_IVAR'.format(linename)] = boxflux_ivar # * u.second**2*u.cm**4/u.erg**2 + + result['{}_CHI2'.format(linename)] = np.sum(emlineivar[lineindx] * (emlineflux[lineindx] - finalmodel[lineindx])**2) # keep track of sigma and z but only using XX-sigma lines linesnr = result['{}_AMP'.format(linename)] * np.sqrt(result['{}_AMP_IVAR'.format(linename)]) @@ -1144,17 +1193,17 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, narrow_redshifts.append(linez) # next, get the continuum, the inverse variance in the line-amplitude, and the EW - indxlo = np.where((emlinewave > (linezwave - 10*linesigma * linezwave / C_LIGHT)) * - (emlinewave < (linezwave - 3.*linesigma * linezwave / C_LIGHT)) * + indxlo = np.where((emlinewave > (linezwave - 10*linesigma_ang_window)) * + (emlinewave < (linezwave - 3.*linesigma_ang_window)) * (oemlineivar > 0))[0] #(finalmodel == 0))[0] - indxhi = np.where((emlinewave < (linezwave + 10*linesigma * linezwave / C_LIGHT)) * - (emlinewave > (linezwave + 3.*linesigma * linezwave / C_LIGHT)) * + indxhi = np.where((emlinewave < (linezwave + 10*linesigma_ang_window)) * + (emlinewave > (linezwave + 3.*linesigma_ang_window)) * (oemlineivar > 0))[0] #(finalmodel == 0))[0] indx = np.hstack((indxlo, indxhi)) - - if len(indx) >= 3: # require at least XX pixels to get the continuum level + + if len(indx) >= nminpix: # require at least XX pixels to get the continuum level #_, cmed, csig = sigma_clipped_stats(specflux_nolines[indx], sigma=3.0) clipflux, _, _ = sigmaclip(specflux_nolines[indx], low=3, high=3) # corner case: if a portion of a camera is masked @@ -1174,6 +1223,7 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, if result['{}_CONT'.format(linename)] != 0.0 and result['{}_CONT_IVAR'.format(linename)] != 0.0: lineflux = result['{}_FLUX'.format(linename)] + #linefluxivar = result['{}_BOXFLUX_IVAR'.format(linename)] linefluxivar = result['{}_FLUX_IVAR'.format(linename)] if lineflux > 0 and linefluxivar > 0: # add the uncertainties in flux and the continuum in quadrature @@ -1192,7 +1242,7 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, result['{}_EW_LIMIT'.format(linename)] = ewlimit if 'debug' in log.name: - for col in ('VSHIFT', 'SIGMA', 'AMP', 'AMP_IVAR', 'CHI2', 'NPIX'): + for col in ('VSHIFT', 'SIGMA', 'MODELAMP', 'AMP', 'AMP_IVAR', 'CHI2', 'NPIX'): log.debug('{} {}: {:.4f}'.format(linename, col, result['{}_{}'.format(linename, col)])) for col in ('FLUX', 'BOXFLUX', 'FLUX_IVAR', 'BOXFLUX_IVAR', 'CONT', 'CONT_IVAR', 'EW', 'EW_IVAR', 'FLUX_LIMIT', 'EW_LIMIT'): log.debug('{} {}: {:.4f}'.format(linename, col, result['{}_{}'.format(linename, col)])) @@ -1228,24 +1278,28 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, # pdb.set_trace() # Clean up the doublets whose amplitudes were tied in the fitting since - # they may have been zeroed out in the clean-up, above. - if result['OIII_5007_AMP'] == 0.0 and result['OIII_5007_NPIX'] > 0: + # they may have been zeroed out in the clean-up, above. This should be + # smarter. + if result['OIII_5007_MODELAMP'] == 0.0 and result['OIII_5007_NPIX'] > 0: + result['OIII_4959_MODELAMP'] = 0.0 result['OIII_4959_AMP'] = 0.0 result['OIII_4959_FLUX'] = 0.0 result['OIII_4959_EW'] = 0.0 - if result['NII_6584_AMP'] == 0.0 and result['NII_6584_NPIX'] > 0: + if result['NII_6584_MODELAMP'] == 0.0 and result['NII_6584_NPIX'] > 0: + result['NII_6548_MODELAMP'] = 0.0 result['NII_6548_AMP'] = 0.0 result['NII_6548_FLUX'] = 0.0 result['NII_6548_EW'] = 0.0 - if result['OII_7320_AMP'] == 0.0 and result['OII_7320_NPIX'] > 0: + if result['OII_7320_MODELAMP'] == 0.0 and result['OII_7320_NPIX'] > 0: + result['OII_7330_MODELAMP'] = 0.0 result['OII_7330_AMP'] = 0.0 result['OII_7330_FLUX'] = 0.0 result['OII_7330_EW'] = 0.0 - if result['MGII_2796_AMP'] == 0.0 and result['MGII_2803_AMP'] == 0.0: + if result['MGII_2796_MODELAMP'] == 0.0 and result['MGII_2803_MODELAMP'] == 0.0: result['MGII_DOUBLET_RATIO'] = 0.0 - if result['OII_3726_AMP'] == 0.0 and result['OII_3729_AMP'] == 0.0: + if result['OII_3726_MODELAMP'] == 0.0 and result['OII_3729_MODELAMP'] == 0.0: result['OII_DOUBLET_RATIO'] = 0.0 - if result['SII_6716_AMP'] == 0.0 and result['SII_6731_AMP'] == 0.0: + if result['SII_6716_MODELAMP'] == 0.0 and result['SII_6731_MODELAMP'] == 0.0: result['SII_DOUBLET_RATIO'] = 0.0 if 'debug' in log.name: @@ -1652,7 +1706,7 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', amp = fastspec['{}_AMP'.format(linename)] if amp != 0: desiemlines_oneline1 = build_emline_model( - self.log10wave, redshift, np.array([amp]), + self.dlog10wave, redshift, np.array([amp]), np.array([fastspec['{}_VSHIFT'.format(linename)]]), np.array([fastspec['{}_SIGMA'.format(linename)]]), np.array([oneline['restwave']]), data['wave'], data['res']) @@ -2068,9 +2122,6 @@ def major_formatter(x, pos): emlinesigma = emlinesigma[good] emlinemodel = emlinemodel[good] - #if icam == 0: - # import matplotlib.pyplot as plt ; plt.clf() ; plt.plot(emlinewave, emlineflux) ; plt.plot(emlinewave, emlinemodel) ; plt.xlim(4180, 4210) ; plt.ylim(-15, 17) ; plt.savefig('desi-users/ioannis/tmp/junkg.png') - emlinemodel_oneline = [] for desiemlines_oneline1 in desiemlines_oneline: emlinemodel_oneline.append(desiemlines_oneline1[icam][good]) @@ -2347,7 +2398,7 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum t0 = time.time() initfit = EMFit.optimize(initial_linemodel_nobroad, emlinewave, emlineflux, weights, redshift, resolution_matrix, camerapix, - log=log, debug=False) + log=log, debug=False, get_finalamp=False) initmodel = EMFit.bestfit(initfit, redshift, emlinewave, resolution_matrix, camerapix) initchi2 = EMFit.chi2(initfit, emlinewave, emlineflux, emlineivar, initmodel) nfree = np.sum((initfit['fixed'] == False) * (initfit['tiedtoparam'] == -1)) @@ -2391,13 +2442,13 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum t0 = time.time() broadfit = EMFit.optimize(initial_linemodel, emlinewave, emlineflux, weights, - redshift, resolution_matrix, camerapix, log=log, debug=False) + redshift, resolution_matrix, camerapix, log=log, + debug=False, get_finalamp=True) broadmodel = EMFit.bestfit(broadfit, redshift, emlinewave, resolution_matrix, camerapix) broadchi2 = EMFit.chi2(broadfit, emlinewave, emlineflux, emlineivar, broadmodel) nfree = np.sum((broadfit['fixed'] == False) * (broadfit['tiedtoparam'] == -1)) log.info('Second (broad) line-fitting with {} free parameters took {:.2f} seconds [niter={}, rchi2={:.4f}].'.format( nfree, time.time()-t0, broadfit.meta['nfev'], broadchi2)) - linechi2_broad, linechi2_init = broadchi2, initchi2 log.info('Chi2 with broad lines = {:.5f} and without broad lines = {:.5f} [chi2_narrow-chi2_broad={:.5f}]'.format( @@ -2415,7 +2466,8 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum sigdrop1 = (broadfit[Habroad]['value'] <= broadfit[broadfit['param_name'] == 'halpha_sigma']['value'])[0] sigdrop2 = broadfit[Habroad]['value'][0] < EMFit.minsigma_balmer_broad - ampsnr = broadfit[Bbroad]['value'].data * np.sqrt(broadfit[Bbroad]['civar'].data) + ampsnr = broadfit[Bbroad]['obsvalue'].data * np.sqrt(broadfit[Bbroad]['civar'].data) + #ampdrop = np.any(ampsnr[-1:] < EMFit.minsnr_balmer_broad) ampdrop = np.any(ampsnr[-2:] < EMFit.minsnr_balmer_broad) #W = (initfit['fixed'] == False) * (initfit['tiedtoparam']==-1) @@ -2527,7 +2579,7 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum t0 = time.time() finalfit = EMFit.optimize(linemodel, emlinewave, emlineflux, weights, redshift, resolution_matrix, camerapix, - log=log, debug=False) + log=log, debug=False, get_finalamp=True) finalmodel = EMFit.bestfit(finalfit, redshift, emlinewave, resolution_matrix, camerapix) finalchi2 = EMFit.chi2(finalfit, emlinewave, emlineflux, emlineivar, finalmodel) nfree = np.sum((finalfit['fixed'] == False) * (finalfit['tiedtoparam'] == -1)) diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index 4eed519e..d77d076e 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -40,7 +40,8 @@ def _assign_units_to_columns(fastfit, metadata, Spec, templates, fastphot, stack def fastspec_one(iobj, data, out, meta, templates, log=None, minspecwave=3500., maxspecwave=9900., broadlinefit=True, - fastphot=False, stackfit=False, nophoto=False, percamera_models=False): + fastphot=False, stackfit=False, nophoto=False, + percamera_models=False, debug_plots=False): """Multiprocessing wrapper to run :func:`fastspec` on a single object. """ @@ -60,7 +61,7 @@ def fastspec_one(iobj, data, out, meta, templates, log=None, continuummodel, smooth_continuum = continuum_specfit(data, out, templatecache, fastphot=fastphot, nophoto=nophoto, - log=log) + debug_plots=debug_plots, log=log) # Optionally fit the emission-line spectrum. if fastphot: @@ -125,6 +126,7 @@ def parse(options=None, log=None): parser.add_argument('--mapdir', type=str, default=None, help='Optional directory name for the dust maps.') parser.add_argument('--dr9dir', type=str, default=None, help='Optional directory name for the DR9 photometry.') parser.add_argument('--specproddir', type=str, default=None, help='Optional directory name for the spectroscopic production.') + parser.add_argument('--debug-plots', action='store_true', help='Generate a variety of debugging plots (written to $PWD).') parser.add_argument('--verbose', action='store_true', help='Be verbose (for debugging purposes).') if log is None: @@ -225,7 +227,7 @@ def fastspec(fastphot=False, stackfit=False, args=None, comm=None, verbose=False t0 = time.time() fitargs = [(iobj, data[iobj], out[iobj], meta[iobj], templates, log, minspecwave, maxspecwave, args.broadlinefit, fastphot, stackfit, - args.nophoto, args.percamera_models) + args.nophoto, args.percamera_models, args.debug_plots) for iobj in np.arange(Spec.ntargets)] if args.mp > 1: import multiprocessing @@ -640,10 +642,10 @@ def major_formatter(x, pos): for oneline in EMFit.linetable[inrange]: # for all lines in range linename = oneline['name'].upper() - amp = fastspec['{}_AMP'.format(linename)] + amp = fastspec['{}_MODELAMP'.format(linename)] if amp != 0: desiemlines_oneline1 = build_emline_model( - EMFit.log10wave, redshift, np.array([amp]), + EMFit.dlog10wave, redshift, np.array([amp]), np.array([fastspec['{}_VSHIFT'.format(linename)]]), np.array([fastspec['{}_SIGMA'.format(linename)]]), np.array([oneline['restwave']]), data['wave'], data['res']) @@ -1016,11 +1018,13 @@ def major_formatter(x, pos): plotsig_default_balmer = 500.0 # [km/s] plotsig_default_broad = 2000.0 # [km/s] - meanwaves, deltawaves, sigmas, linenames = [], [], [], [] + minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames = [], [], [], [], [], [] for plotgroup in set(linetable['plotgroup']): I = np.where(plotgroup == linetable['plotgroup'])[0] linename = linetable['nicename'][I[0]].replace('-', ' ') linenames.append(linename) + minwaves.append(np.min(linetable['restwave'][I])) + maxwaves.append(np.max(linetable['restwave'][I])) meanwaves.append(np.mean(linetable['restwave'][I])) deltawaves.append((np.max(linetable['restwave'][I]) - np.min(linetable['restwave'][I])) / 2) @@ -1053,6 +1057,8 @@ def major_formatter(x, pos): ax = [] else: srt = np.argsort(meanwaves) + minwaves = np.hstack(minwaves)[srt] + maxwaves = np.hstack(maxwaves)[srt] meanwaves = np.hstack(meanwaves)[srt] deltawaves = np.hstack(deltawaves)[srt] sigmas = np.hstack(sigmas)[srt] @@ -1083,7 +1089,8 @@ def major_formatter(x, pos): else: ax, irow, colshift = [], 4, 5 # skip the gap row - for iax, (meanwave, deltawave, sig, linename) in enumerate(zip(meanwaves, deltawaves, sigmas, linenames)): + for iax, (minwave, maxwave, meanwave, deltawave, sig, linename) in enumerate( + zip(minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames)): icol = iax % nlinecols icol += colshift if iax > 0 and iax % nlinecols == 0: @@ -1093,9 +1100,10 @@ def major_formatter(x, pos): xx = fig.add_subplot(gs[irow, icol]) ax.append(xx) - wmin = (meanwave - deltawave) * (1+redshift) - 6 * sig * meanwave * (1+redshift) / C_LIGHT - wmax = (meanwave + deltawave) * (1+redshift) + 6 * sig * meanwave * (1+redshift) / C_LIGHT - #print(linename, wmin, wmax) + wmin = (minwave - deltawave) * (1+redshift) - 6 * sig * minwave * (1+redshift) / C_LIGHT + wmax = (maxwave + deltawave) * (1+redshift) + 6 * sig * maxwave * (1+redshift) / C_LIGHT + #wmin = (meanwave - deltawave) * (1+redshift) - 6 * sig * meanwave * (1+redshift) / C_LIGHT + #wmax = (meanwave + deltawave) * (1+redshift) + 6 * sig * meanwave * (1+redshift) / C_LIGHT # iterate over cameras for icam in np.arange(len(data['cameras'])): # iterate over cameras @@ -1109,9 +1117,6 @@ def major_formatter(x, pos): emlinesigma = emlinesigma[good] emlinemodel = emlinemodel[good] - #if icam == 0: - # import matplotlib.pyplot as plt ; plt.clf() ; plt.plot(emlinewave, emlineflux) ; plt.plot(emlinewave, emlinemodel) ; plt.xlim(4180, 4210) ; plt.ylim(-15, 17) ; plt.savefig('desi-users/ioannis/tmp/junkg.png') - emlinemodel_oneline = [] for desiemlines_oneline1 in desiemlines_oneline: emlinemodel_oneline.append(desiemlines_oneline1[icam][good]) diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index 829656be..ab52a18a 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -101,22 +101,25 @@ def _unpack_one_spectrum(args): """Multiprocessing wrapper.""" return unpack_one_spectrum(*args) -def unpack_one_spectrum(iobj, specdata, meta, ebv, Filters, fastphot, synthphot, log): +def unpack_one_spectrum(iobj, specdata, meta, ebv, fastphot, synthphot, log): """Unpack the data for a single object and correct for Galactic extinction. Also flag pixels which may be affected by emission lines. """ + from fastspecfit.continuum import ContinuumTools from desiutil.dust import mwdust_transmission, dust_transmission log.info('Pre-processing object {} [targetid {} z={:.6f}].'.format( iobj, meta['TARGETID'], meta['Z'])) + CTools = ContinuumTools() + if specdata['photsys'] == 'S': - filters = Filters.decam - allfilters = Filters.decamwise + filters = CTools.decam + allfilters = CTools.decamwise else: - filters = Filters.bassmzls - allfilters = Filters.bassmzlswise + filters = CTools.bassmzls + allfilters = CTools.bassmzlswise RV = 3.1 @@ -127,16 +130,16 @@ def unpack_one_spectrum(iobj, specdata, meta, ebv, Filters, fastphot, synthphot, # self-consistent with how we correct the photometry for dust. meta['EBV'] = ebv if specdata['photsys'] != '': - mw_transmission_flux = np.array([mwdust_transmission(ebv, band, specdata['photsys'], match_legacy_surveys=False) for band in Filters.bands]) - for band, mwdust in zip(Filters.bands, mw_transmission_flux): + mw_transmission_flux = np.array([mwdust_transmission(ebv, band, specdata['photsys'], match_legacy_surveys=False) for band in CTools.bands]) + for band, mwdust in zip(CTools.bands, mw_transmission_flux): meta['MW_TRANSMISSION_{}'.format(band.upper())] = mwdust else: - #mw_transmission_fiberflux = np.ones(len(Filters.bands)) + #mw_transmission_fiberflux = np.ones(len(CTools.bands)) mw_transmission_flux = 10**(-0.4 * ebv * RV * ext_odonnell(allfilters.effective_wavelengths.value, Rv=RV)) - maggies = np.zeros(len(Filters.bands)) - ivarmaggies = np.zeros(len(Filters.bands)) - for iband, band in enumerate(Filters.bands): + maggies = np.zeros(len(CTools.bands)) + ivarmaggies = np.zeros(len(CTools.bands)) + for iband, band in enumerate(CTools.bands): maggies[iband] = meta['FLUX_{}'.format(band.upper())] / mw_transmission_flux[iband] ivarmaggies[iband] = meta['FLUX_IVAR_{}'.format(band.upper())] * mw_transmission_flux[iband]**2 @@ -145,30 +148,30 @@ def unpack_one_spectrum(iobj, specdata, meta, ebv, Filters, fastphot, synthphot, log.critical(errmsg) raise ValueError(errmsg) - specdata['phot'] = Filters.parse_photometry( - Filters.bands, maggies=maggies, ivarmaggies=ivarmaggies, nanomaggies=True, + specdata['phot'] = CTools.parse_photometry( + CTools.bands, maggies=maggies, ivarmaggies=ivarmaggies, nanomaggies=True, lambda_eff=allfilters.effective_wavelengths.value, - min_uncertainty=Filters.min_uncertainty, log=log) + min_uncertainty=CTools.min_uncertainty, log=log) # fiber fluxes if specdata['photsys'] != '': - mw_transmission_fiberflux = np.array([mwdust_transmission(ebv, band, specdata['photsys']) for band in Filters.fiber_bands]) + mw_transmission_fiberflux = np.array([mwdust_transmission(ebv, band, specdata['photsys']) for band in CTools.fiber_bands]) else: - #mw_transmission_fiberflux = np.ones(len(Filters.fiber_bands)) + #mw_transmission_fiberflux = np.ones(len(CTools.fiber_bands)) mw_transmission_fiberflux = 10**(-0.4 * ebv * RV * ext_odonnell(filters.effective_wavelengths.value, Rv=RV)) - fibermaggies = np.zeros(len(Filters.fiber_bands)) - fibertotmaggies = np.zeros(len(Filters.fiber_bands)) - #ivarfibermaggies = np.zeros(len(Filters.fiber_bands)) - for iband, band in enumerate(Filters.fiber_bands): + fibermaggies = np.zeros(len(CTools.fiber_bands)) + fibertotmaggies = np.zeros(len(CTools.fiber_bands)) + #ivarfibermaggies = np.zeros(len(CTools.fiber_bands)) + for iband, band in enumerate(CTools.fiber_bands): fibermaggies[iband] = meta['FIBERFLUX_{}'.format(band.upper())] / mw_transmission_fiberflux[iband] fibertotmaggies[iband] = meta['FIBERTOTFLUX_{}'.format(band.upper())] / mw_transmission_fiberflux[iband] #ivarfibermaggies[iband] = meta['FIBERTOTFLUX_IVAR_{}'.format(band.upper())] * mw_transmission_fiberflux[iband]**2 - specdata['fiberphot'] = Filters.parse_photometry(Filters.fiber_bands, + specdata['fiberphot'] = CTools.parse_photometry(CTools.fiber_bands, maggies=fibermaggies, nanomaggies=True, lambda_eff=filters.effective_wavelengths.value, log=log) - specdata['fibertotphot'] = Filters.parse_photometry(Filters.fiber_bands, + specdata['fibertotphot'] = CTools.parse_photometry(CTools.fiber_bands, maggies=fibertotmaggies, nanomaggies=True, lambda_eff=filters.effective_wavelengths.value, log=log) @@ -233,9 +236,9 @@ def unpack_one_spectrum(iobj, specdata, meta, ebv, Filters, fastphot, synthphot, specdata['camerapix'][icam, :] = [np.sum(npixpercam[:icam+1]), np.sum(npixpercam[:icam+2])] # coadded spectrum - coadd_linemask_dict = Filters.build_linemask(specdata['coadd_wave'], specdata['coadd_flux'], - specdata['coadd_ivar'], redshift=specdata['zredrock'], - linetable=Filters.linetable) + coadd_linemask_dict = CTools.build_linemask(specdata['coadd_wave'], specdata['coadd_flux'], + specdata['coadd_ivar'], redshift=specdata['zredrock'], + linetable=CTools.linetable) specdata['coadd_linename'] = coadd_linemask_dict['linename'] specdata['coadd_linepix'] = [np.where(lpix)[0] for lpix in coadd_linemask_dict['linepix']] specdata['coadd_contpix'] = [np.where(cpix)[0] for cpix in coadd_linemask_dict['contpix']] @@ -293,11 +296,11 @@ def unpack_one_spectrum(iobj, specdata, meta, ebv, Filters, fastphot, synthphot, #synthvarmaggies = filters.get_ab_maggies(1e-17**2 * padvar, padwave) #synthivarmaggies = 1 / synthvarmaggies.as_array().view('f8')[:3] # keep just grz - #specdata['synthphot'] = Filters.parse_photometry(Filters.bands, + #specdata['synthphot'] = CTools.parse_photometry(CTools.bands, # maggies=synthmaggies, lambda_eff=lambda_eff[:3], # ivarmaggies=synthivarmaggies, nanomaggies=False, log=log) - specdata['synthphot'] = Filters.parse_photometry(Filters.synth_bands, + specdata['synthphot'] = CTools.parse_photometry(CTools.synth_bands, maggies=synthmaggies, nanomaggies=False, lambda_eff=filters.effective_wavelengths.value, log=log) @@ -307,29 +310,33 @@ def _unpack_one_stacked_spectrum(args): """Multiprocessing wrapper.""" return unpack_one_stacked_spectrum(*args) -def unpack_one_stacked_spectrum(iobj, specdata, meta, Filters, synthphot, log): +def unpack_one_stacked_spectrum(iobj, specdata, meta, synthphot, log): """Unpack the data for a single stacked spectrum. Also flag pixels which may be affected by emission lines. """ + from fastspecfit.continuum import ContinuumTools + + CTools = ContinuumTools() + log.info('Pre-processing object {} [stackid {} z={:.6f}].'.format( iobj, meta['STACKID'], meta['Z'])) if specdata['photsys'] == 'S': - filters = Filters.decam - allfilters = Filters.decamwise + filters = CTools.decam + allfilters = CTools.decamwise else: - filters = Filters.bassmzls - allfilters = Filters.bassmzlswise + filters = CTools.bassmzls + allfilters = CTools.bassmzlswise # Dummy imaging photometry. - maggies = np.zeros(len(Filters.bands)) - ivarmaggies = np.zeros(len(Filters.bands)) + maggies = np.zeros(len(CTools.bands)) + ivarmaggies = np.zeros(len(CTools.bands)) - specdata['phot'] = Filters.parse_photometry( - Filters.bands, maggies=maggies, ivarmaggies=ivarmaggies, nanomaggies=True, + specdata['phot'] = CTools.parse_photometry( + CTools.bands, maggies=maggies, ivarmaggies=ivarmaggies, nanomaggies=True, lambda_eff=allfilters.effective_wavelengths.value, - min_uncertainty=Filters.min_uncertainty, log=log) + min_uncertainty=CTools.min_uncertainty, log=log) #specdata['fiberphot'] = specdata['phot'] #specdata['fibertotphot'] = specdata['phot'] @@ -391,9 +398,9 @@ def unpack_one_stacked_spectrum(iobj, specdata, meta, Filters, synthphot, log): del specdata[key] # coadded spectrum - coadd_linemask_dict = Filters.build_linemask(specdata['coadd_wave'], specdata['coadd_flux'], - specdata['coadd_ivar'], redshift=specdata['zredrock'], - linetable=Filters.linetable) + coadd_linemask_dict = CTools.build_linemask(specdata['coadd_wave'], specdata['coadd_flux'], + specdata['coadd_ivar'], redshift=specdata['zredrock'], + linetable=CTools.linetable) specdata['coadd_linename'] = coadd_linemask_dict['linename'] specdata['coadd_linepix'] = [np.where(lpix)[0] for lpix in coadd_linemask_dict['linepix']] specdata['coadd_contpix'] = [np.where(cpix)[0] for cpix in coadd_linemask_dict['contpix']] @@ -437,7 +444,7 @@ def unpack_one_stacked_spectrum(iobj, specdata, meta, Filters, synthphot, log): synthmaggies = filters.get_ab_maggies(padflux / FLUXNORM, padwave) synthmaggies = synthmaggies.as_array().view('f8') - specdata['synthphot'] = Filters.parse_photometry(Filters.synth_bands, + specdata['synthphot'] = CTools.parse_photometry(CTools.synth_bands, maggies=synthmaggies, nanomaggies=False, lambda_eff=filters.effective_wavelengths.value, log=log) @@ -1087,9 +1094,7 @@ def read_and_unpack(self, fastphot=False, synthphot=True, mp=1): from desispec.coaddition import coadd_cameras from desispec.io import read_spectra from desiutil.dust import SFDMap - from fastspecfit.continuum import ContinuumTools - CTools = ContinuumTools() SFD = SFDMap(scaling=1.0, mapdir=self.mapdir) alldata = [] @@ -1115,7 +1120,7 @@ def read_and_unpack(self, fastphot=False, synthphot=True, mp=1): 'photsys': meta['PHOTSYS'][iobj], 'dluminosity': dlum[iobj], 'dmodulus': dmod[iobj], 'tuniv': tuniv[iobj], } - unpackargs.append((iobj, specdata, meta[iobj], ebv[iobj], CTools, True, False, log)) + unpackargs.append((iobj, specdata, meta[iobj], ebv[iobj], True, False, log)) else: from desispec.resolution import Resolution @@ -1147,7 +1152,7 @@ def read_and_unpack(self, fastphot=False, synthphot=True, mp=1): 'coadd_ivar': coadd_spec.ivar[coadd_cameras][iobj, :], 'coadd_res': Resolution(coadd_spec.resolution_data[coadd_cameras][iobj, :]), } - unpackargs.append((iobj, specdata, meta[iobj], ebv[iobj], CTools, fastphot, synthphot, log)) + unpackargs.append((iobj, specdata, meta[iobj], ebv[iobj], fastphot, synthphot, log)) if mp > 1: import multiprocessing @@ -1342,8 +1347,6 @@ def read_stacked(self, stackfiles, firsttarget=0, ntargets=None, return # Now read the data as in self.read_and_unpack (for unstacked spectra). - CTools = ContinuumTools() - alldata = [] for ispec, (stackfile, meta) in enumerate(zip(self.stackfiles, self.meta)): nobj = len(meta) @@ -1397,7 +1400,7 @@ def read_stacked(self, stackfiles, firsttarget=0, ntargets=None, 'coadd_ivar': specdata['ivar0'][0], 'coadd_res': specdata['res0'][0], }) - unpackargs.append((iobj, specdata, meta[iobj], CTools, synthphot, log)) + unpackargs.append((iobj, specdata, meta[iobj], synthphot, log)) if mp > 1: import multiprocessing @@ -1661,6 +1664,8 @@ def init_fastspec_output(input_meta, specprod, templates=None, ncoeff=None, for line in linetable['name']: line = line.upper() + out.add_column(Column(name='{}_MODELAMP'.format(line), length=nobj, dtype='f4', + unit=10**(-17)*u.erg/(u.second*u.cm**2*u.Angstrom))) out.add_column(Column(name='{}_AMP'.format(line), length=nobj, dtype='f4', unit=10**(-17)*u.erg/(u.second*u.cm**2*u.Angstrom))) out.add_column(Column(name='{}_AMP_IVAR'.format(line), length=nobj, dtype='f4', @@ -1737,7 +1742,7 @@ def read_fastspecfit(fastfitfile, rows=None, columns=None, read_models=False): # Add specprod to the metadata table so that we can stack across # productions (e.g., Fuji+Guadalupe). - hdr = fitsio.read_header(fastfitfile, ext='PRIMARY') + hdr = fitsio.read_header(fastfitfile, ext=0) if 'SPECPROD' in hdr: specprod = hdr['SPECPROD'] meta['SPECPROD'] = specprod diff --git a/py/fastspecfit/util.py b/py/fastspecfit/util.py index ee08a6f3..558d6d31 100644 --- a/py/fastspecfit/util.py +++ b/py/fastspecfit/util.py @@ -231,7 +231,7 @@ def find_minima(x): return ii[jj] -def minfit(x, y): +def minfit(x, y, return_coeff=False): """Fits y = y0 + ((x-x0)/xerr)**2 See redrock.zwarning.ZWarningMask.BAD_MINFIT for zwarn failure flags @@ -244,17 +244,27 @@ def minfit(x, y): (tuple): (x0, xerr, y0, zwarn) where zwarn=0 is good fit. """ + a, b, c = 0., 0., 0. if len(x) < 3: - return (-1,-1,-1,ZWarningMask.BAD_MINFIT) + if return_coeff: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c)) + else: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT) try: #- y = a x^2 + b x + c a,b,c = np.polyfit(x,y,2) except np.linalg.LinAlgError: - return (-1,-1,-1,ZWarningMask.BAD_MINFIT) + if return_coeff: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c)) + else: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT) if a == 0.0: - return (-1,-1,-1,ZWarningMask.BAD_MINFIT) + if return_coeff: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c)) + else: + return (-1,-1,-1,ZWarningMask.BAD_MINFIT) #- recast as y = y0 + ((x-x0)/xerr)^2 x0 = -b / (2*a) @@ -272,7 +282,10 @@ def minfit(x, y): xerr = 1 / np.sqrt(-a) zwarn |= ZWarningMask.BAD_MINFIT - return (x0, xerr, y0, zwarn) + if return_coeff: + return (x0, xerr, y0, zwarn,(a,b,c)) + else: + return (x0, xerr, y0, zwarn) class TabulatedDESI(object): """ @@ -292,8 +305,8 @@ class TabulatedDESI(object): """ def __init__(self): - from pkg_resources import resource_filename - cosmofile = resource_filename('fastspecfit', 'data/desi_fiducial_cosmology.dat') + from importlib import resources + cosmofile = resources.files('fastspecfit').joinpath('data/desi_fiducial_cosmology.dat') self._z, self._efunc, self._comoving_radial_distance = np.loadtxt(cosmofile, comments='#', usecols=None, unpack=True) diff --git a/sandbox/nnls.py b/sandbox/nnls.py new file mode 100644 index 00000000..51e8edad --- /dev/null +++ b/sandbox/nnls.py @@ -0,0 +1,55 @@ +#import numba +#@numba.jit(nopython=True) +def nnls(A, b, maxiter=None, eps=1e-7, v=False): + # https://en.wikipedia.org/wiki/Non-negative_least_squares#Algorithms + # not sure what eps should be set to + m, n = A.shape + if b.shape[0] != m: + raise ValueError()#f"Shape mismatch: {A.shape} {b.shape}. Expected: (m, n) (m) ") + if maxiter is None: + # this is what scipy does when maxiter is not specified + maxiter = 3 * n + # init + P = np.zeros(n).astype(bool) + R = np.ones(n).astype(bool) + x = np.zeros(n) + s = np.zeros(n) + # compute these once ahead of time + ATA = A.T.dot(A) + ATb = A.T.dot(b) + # main loop + w = ATb - ATA.dot(x) + j = np.argmax(R * w) + c = 0 # iteration count + # while R != {} and max(w[R]) > eps + while np.any(R) and w[j] > eps: + #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}") + # add j to P, remove j from R + P[j], R[j] = True, False + s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P]) + s[R] = 0 + d = 0 # inner loop iteration count, for debugging + #if v: print(f"{c=}", f"{P=}\n {R=}\n {s=}") + # while P != {} and min(s[P]) < eps + # make sure P is not empty before checking min s[P] + while np.any(P) and np.min(s[P]) < eps: + i = P & (s < eps) + #if v: print(f" {d=}", f" {P=}\n {i=}\n {s=}\n {x=}") + a = np.min(x[i] / (x[i] - s[i])) + x = x + a * (s - x) + j = P & (x < eps) + R[j], P[j] = True, False + s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P]) + s[R] = 0 + d += 1 + x[:] = s + # w = A.T.dot(b - A.dot(x)) + w = ATb - ATA.dot(x) + j = np.argmax(R * w) + #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}") + c += 1 + if c >= maxiter: + break + res = np.linalg.norm(A.dot(x) - b) + return x, res +