diff --git a/.gitignore b/.gitignore index 19b9258..41d3a46 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,14 @@ build/ dist/ multiprofit.egg-info/ -*/__pycache__/ +__pycache__/ +.coverage +.sconsign.dblite +config.log +version.py + +# Generated by pip install -e . +python/lsst_multiprofit.egg-info/ # IDE folders and files .cache/ @@ -11,6 +18,9 @@ multiprofit.egg-info/ .vscode/ compile_commands.json +# test outputs +tests/.tests + # pytest plugins .hypothesis/ prof/ diff --git a/SConstruct b/SConstruct new file mode 100644 index 0000000..1ad1023 --- /dev/null +++ b/SConstruct @@ -0,0 +1,4 @@ +# -*- python -*- +from lsst.sconsUtils import scripts +# Python-only package +scripts.BasicSConstruct("multiprofit", disableCc=True, noCfgFile=True) diff --git a/examples/fithsc.py b/examples/fithsc.py index 52161e7..9ae9c2d 100644 --- a/examples/fithsc.py +++ b/examples/fithsc.py @@ -18,9 +18,9 @@ import gauss2d.fit as g2f import matplotlib as mpl import matplotlib.pyplot as plt -from multiprofit.componentconfig import SersicConfig, SersicIndexConfig -from multiprofit.fit_psf import CatalogExposurePsfABC, CatalogPsfFitterConfig, CatalogPsfFitter -from multiprofit.fit_source import CatalogExposureSourcesABC, CatalogSourceFitterABC, CatalogSourceFitterConfig +from lsst.multiprofit.componentconfig import SersicConfig, SersicIndexConfig +from lsst.multiprofit.fit_psf import CatalogExposurePsfABC, CatalogPsfFitterConfig, CatalogPsfFitter +from lsst.multiprofit.fit_source import CatalogExposureSourcesABC, CatalogSourceFitterABC, CatalogSourceFitterConfig import numpy as np from typing import Any, Iterable, Mapping diff --git a/examples/test_gaussians.py b/examples/test_gaussians.py index f3c6923..1a48893 100644 --- a/examples/test_gaussians.py +++ b/examples/test_gaussians.py @@ -6,7 +6,7 @@ nbenchmark=100, do_like=True, do_residual=True, do_grad=True, do_jac=True, do_meas_modelfit=False, nsub=4, ) -for x in test: +for x in test: print(f"re={x['reff']:.3f} q={x['axrat']:.2f}" f" ang={x['ang']:2.1f} { x['string']}") print(f'Test complete in {timer() - start:.2f}s') diff --git a/examples/test_utils.py b/examples/test_utils.py index 59ad519..09b8cc2 100644 --- a/examples/test_utils.py +++ b/examples/test_utils.py @@ -22,7 +22,8 @@ from importlib.util import find_spec import numpy as np import matplotlib.pyplot as plt -from multiprofit.fitutils import get_model +# TODO: Replace with e.g. fit_source.get_model? +from lsst.multiprofit.fitutils import get_model import gauss2d as g2 import timeit @@ -181,7 +182,7 @@ def get_setup(xdim=15, ydim=15, r_major=1, axrat=1, angle=0, nsub=1, f'xdim={xdim}', f'ydim={ydim}', f'centroid = g2.Centroid({xdim}/2, {ydim}/2)', - f'kernel = g2.Gaussian(centroid=g2.Centroid(0, 0),' + 'kernel = g2.Gaussian(centroid=g2.Centroid(0, 0),' 'ellipse=g2.Ellipse(sigma_x=0., sigma_y=0))', 'ellipse = g2.Ellipse(g2.EllipseMajor(' f'r_major={r_major}*g2.M_SIGMA_HWHM, axrat={axrat}, angle={angle}, degrees=True))', @@ -318,7 +319,7 @@ def gaussian_test(xdim=49, ydim=51, reffs=None, angs=None, axrats=None, nbenchma ang_rad = np.deg2rad(ang) for key in ("dev", "exp"): times[f"mmf-{key}"] = np.min(timeit.repeat( - f"msf.evaluate().addToImage(img)", + "msf.evaluate().addToImage(img)", setup=( f"import numpy as np;" f"from lsst.shapelet import RadialProfile;" @@ -360,7 +361,7 @@ def gradient_test(dimx=5, dimy=4, flux=1e4, reff=2, axrat=0.5, ang=0, bg=1e3, reff_psf=0, axrat_psf=0.95, ang_psf=0, n_psfs=1, printout=False, plot=False): if n_psfs > 1: - raise ValueError(f"n_psfs>1 not yet supported") + raise ValueError("n_psfs>1 not yet supported") cen_x, cen_y = dimx/2., dimy/2. # Keep this in units of sigma, not re==FWHM/2 diff --git a/multiprofit/__init__.py b/multiprofit/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/multiprofit/multigaussianapproxprofile.py b/multiprofit/multigaussianapproxprofile.py deleted file mode 100644 index 265fd74..0000000 --- a/multiprofit/multigaussianapproxprofile.py +++ /dev/null @@ -1,1823 +0,0 @@ -# This file is part of multiprofit. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - -import copy -import galsim as gs -import gauss2d -import gauss2d.fit as g2f -import multiprofit.objects as mpfobj -import multiprofit.utils as mpfutil -import numpy as np -import scipy.interpolate as spinterp -import scipy.optimize as spopt -from typing import Dict - -ln10 = np.log(10) - -def _print_weights_c(weights: Dict): - nl = '\n' - for key, value in weights.items(): - norm_strs = [f"{x:.15e}" for x in value[0]] - value[0][-1] = 1.0 - sum(float(x) for x in norm_strs[:-1]) - norm_strs[-1] = f"{value[0][-1]:.15e}" - value_sum = sum(float(x) for x in norm_strs) - if value_sum != 1: - raise ValueError(f'key={key} value_sum-1={value_sum-1} != 0') - string = f'''{{{key}, {{{nl}{f",{nl}".join(f" {{{x}, {y*gauss2d.M_HWHM_SIGMA}}}" - for x, y in zip(norm_strs, value[1]))}{nl}}}}},''' - print(string) - - -# Compute the fraction of the integrated flux within x for a sum of Gaussians -# x is a length in arbitrary units -# Weightsizes is a list of tuples of the weight (flux) and size (r_eff in the units of x) of each gaussian -# Note that gauss2dint expects x/sigma, but size is re, so we pass x/re*re/sigma = x/sigma -# 0 > quant > 1 turns it into a function that returns zero -# at the value of x containing a fraction quant of the total flux -# This is so you can use root finding algorithms to find x for a given quant (see below) -def multigauss2dint(x, weightsizes, quant=0): - re_to_sigma = np.sqrt(2.*np.log(2.)) - weight_sum_to_x = 0 - weight_sum = 0 - for weight, size in weightsizes: - weight_sum_to_x += weight*(gauss2dint(x/size*re_to_sigma) if size > 0 else 1) - weight_sum += weight - return weight_sum_to_x/weight_sum - quant - - -# Compute x_quant for a sum of Gaussians, where 0=0 & <=1'.format(quant, quant)) - weight_sum_zero_size = 0 - weight_sum = 0 - for weight, size in weightsizes: - if not (size > 0): - weight_sum_zero_size += weight - weight_sum += weight - if weight_sum_zero_size/weight_sum >= quant: - return 0 - return spopt.brentq(multigauss2dint, a=xmin, b=xmax, args=(weightsizes, quant)) - - -class MultiGaussianApproximationComponent(mpfobj.EllipticalComponent): - # TODO: Figure out a better way to split functionality between this and EllipticalComponent - # Maybe have a generic EllipticalComponent and SingleEllipticalComponent? - """ - Class for multi-Gaussian profiles with (generalized) ellipse shape(s). - """ - profiles_avail = ["sersic"] - ENGINES = ["galsim", "libprofit"] - mandatory = { - "sersic": ["n_ser"], - } - - # These weights are derived here: - # https://github.com/lsst-dm/modelling_research/blob/master/jupyternotebooks/multigaussian_sersic1d.ipynb - weights = { - "sersic": { - 4: { - 0.5: ( - mpfutil.normalize( - np.array([1., 0., 0., 0.])), - np.array([1.0000000000e+00, 8.6800000000e-01, 5.5600000000e-01, 2.7200000000e-01]), - ), - 0.501: ( - mpfutil.normalize( - np.array([9.8945700264e-01, 9.9340700388e-03, 5.9654256544e-04, 1.2384756796e-05])), - np.array([1.0018152162e+00, 8.6829891694e-01, 5.5601124505e-01, 2.7194641925e-01]), - ), - 0.502: ( - mpfutil.normalize( - np.array([9.7915845320e-01, 1.9623440160e-02, 1.1927485937e-03, 2.5358041869e-05])), - np.array([1.0036304324e+00, 8.6859783387e-01, 5.5602249010e-01, 2.7189283849e-01]), - ), - 0.503: ( - mpfutil.normalize( - np.array([9.6909614755e-01, 2.9076336691e-02, 1.7886076942e-03, 3.8908060978e-05])), - np.array([1.0054456487e+00, 8.6889675081e-01, 5.5603373515e-01, 2.7183925774e-01]), - ), - 0.504: ( - mpfutil.normalize( - np.array([9.5926219793e-01, 3.8300684239e-02, 2.3840804072e-03, 5.3037419541e-05])), - np.array([1.0072608649e+00, 8.6919566775e-01, 5.5604498020e-01, 2.7178567698e-01]), - ), - 0.505: ( - mpfutil.normalize( - np.array([9.4964907955e-01, 4.7304034601e-02, 2.9791457766e-03, 6.7740074277e-05])), - np.array([1.0090760811e+00, 8.6949458468e-01, 5.5605622525e-01, 2.7173209623e-01]), - ), - 0.506: ( - mpfutil.normalize( - np.array([9.4024957357e-01, 5.6093640927e-02, 3.5737695079e-03, 8.3015990936e-05])), - np.array([1.0108912973e+00, 8.6979350162e-01, 5.5606747030e-01, 2.7167851547e-01]), - ), - 0.507: ( - mpfutil.normalize( - np.array([9.3105678877e-01, 6.4676421037e-02, 4.1679292284e-03, 9.8860967198e-05])), - np.array([1.0127065136e+00, 8.7009241855e-01, 5.5607871535e-01, 2.7162493472e-01]), - ), - 0.508: ( - mpfutil.normalize( - np.array([9.2206414921e-01, 7.3058966723e-02, 4.7616185406e-03, 1.1526552325e-04])), - np.array([1.0145217298e+00, 8.7039133549e-01, 5.5608996040e-01, 2.7157135396e-01]), - ), - 0.509: ( - mpfutil.normalize( - np.array([9.1326529121e-01, 8.1247682677e-02, 5.3547918794e-03, 1.3223423568e-04])), - np.array([1.0163369460e+00, 8.7069025243e-01, 5.5610120545e-01, 2.7151777321e-01]), - ), - 0.51: ( - mpfutil.normalize( - np.array([9.0465415940e-01, 8.9248645803e-02, 5.9474342283e-03, 1.4976056908e-04])), - np.array([1.0181521622e+00, 8.7098916936e-01, 5.5611245050e-01, 2.7146419246e-01]), - ), - 0.512: ( - mpfutil.normalize( - np.array([8.8797205161e-01, 1.0471043423e-01, 7.1310424559e-03, 1.8647171091e-04])), - np.array([1.0217825947e+00, 8.7158700324e-01, 5.5613494059e-01, 2.7135703095e-01]), - ), - 0.514: ( - mpfutil.normalize( - np.array([8.7197402933e-01, 1.1948833341e-01, 8.3122605259e-03, 2.2537673275e-04])), - np.array([1.0254130271e+00, 8.7218483711e-01, 5.5615743069e-01, 2.7124986944e-01]), - ), - 0.516: ( - mpfutil.normalize( - np.array([8.5661973094e-01, 1.3362289316e-01, 9.4909271779e-03, 2.6644872161e-04])), - np.array([1.0290434596e+00, 8.7278267098e-01, 5.5617992079e-01, 2.7114270793e-01]), - ), - 0.518: ( - mpfutil.normalize( - np.array([8.4187191031e-01, 1.4715152809e-01, 1.0666905708e-02, 3.0965589274e-04])), - np.array([1.0326738920e+00, 8.7338050486e-01, 5.5620241089e-01, 2.7103554642e-01]), - ), - 0.52: ( - mpfutil.normalize( - np.array([8.2769605865e-01, 1.6010893375e-01, 1.1840026965e-02, 3.5498064219e-04])), - np.array([1.0363043244e+00, 8.7397833873e-01, 5.5622490099e-01, 2.7092838491e-01]), - ), - 0.525: ( - mpfutil.normalize( - np.array([7.9455497749e-01, 1.9020815824e-01, 1.4759495284e-02, 4.7736898718e-04])), - np.array([1.0453804056e+00, 8.7547292341e-01, 5.5628112624e-01, 2.7066048114e-01]), - ), - 0.53: ( - mpfutil.normalize( - np.array([7.6435585744e-01, 2.1737336399e-01, 1.7658366124e-02, 6.1241244531e-04])), - np.array([1.0544564867e+00, 8.7696750809e-01, 5.5633735149e-01, 2.7039257736e-01]), - ), - 0.535: ( - mpfutil.normalize( - np.array([7.3673150348e-01, 2.4197382835e-01, 2.0534948388e-02, 7.5971977755e-04])), - np.array([1.0635325678e+00, 8.7846209277e-01, 5.5639357673e-01, 2.7012467359e-01]), - ), - 0.54: ( - mpfutil.normalize( - np.array([7.1137288409e-01, 2.6432043281e-01, 2.3387773784e-02, 9.1890931808e-04])), - np.array([1.0726086489e+00, 8.7995667746e-01, 5.5644980198e-01, 2.6985676982e-01]), - ), - 0.545: ( - mpfutil.normalize( - np.array([6.8801812381e-01, 2.8467668137e-01, 2.6215589953e-02, 1.0896048655e-03])), - np.array([1.0816847300e+00, 8.8145126214e-01, 5.5650602723e-01, 2.6958886605e-01]), - ), - 0.55: ( - mpfutil.normalize( - np.array([6.6644385745e-01, 3.0326737752e-01, 2.9017327056e-02, 1.2714379704e-03])), - np.array([1.0907608111e+00, 8.8294584682e-01, 5.5656225248e-01, 2.6932096228e-01]), - ), - 0.56: ( - mpfutil.normalize( - np.array([6.2789648133e-01, 3.3589734166e-01, 3.4539094688e-02, 1.6670823192e-03])), - np.array([1.1089129733e+00, 8.8593501618e-01, 5.5667470297e-01, 2.6878515473e-01]), - ), - 0.57: ( - mpfutil.normalize( - np.array([5.9448786957e-01, 3.6346156025e-01, 3.9947517585e-02, 2.1030525873e-03])), - np.array([1.1270651355e+00, 8.8892418555e-01, 5.5678715347e-01, 2.6824934719e-01]), - ), - 0.58: ( - mpfutil.normalize( - np.array([5.6527551984e-01, 3.8690886198e-01, 4.5238935286e-02, 2.5766828937e-03])), - np.array([1.1452172978e+00, 8.9191335491e-01, 5.5689960396e-01, 2.6771353964e-01]), - ), - 0.59: ( - mpfutil.normalize( - np.array([5.3953208499e-01, 4.0697128630e-01, 5.0411196453e-02, 3.0854322571e-03])), - np.array([1.1633694600e+00, 8.9490252428e-01, 5.5701205446e-01, 2.6717773209e-01]), - ), - 0.6: ( - mpfutil.normalize( - np.array([5.1668620161e-01, 4.2422360341e-01, 5.5463315784e-02, 3.6268791954e-03])), - np.array([1.1815216222e+00, 8.9789169364e-01, 5.5712450495e-01, 2.6664192455e-01]), - ), - 0.625: ( - mpfutil.normalize( - np.array([4.7758685902e-01, 4.4871014896e-01, 6.8837799255e-02, 4.8651927687e-03])), - np.array([1.2213388674e+00, 9.0522947018e-01, 5.5593716167e-01, 2.6448724702e-01]), - ), - 0.65: ( - mpfutil.normalize( - np.array([4.4790973196e-01, 4.6472790091e-01, 8.1183680488e-02, 6.1786866405e-03])), - np.array([1.2597206861e+00, 9.1227538845e-01, 5.5420949557e-01, 2.6202170433e-01]), - ), - 0.675: ( - mpfutil.normalize( - np.array([4.2468720925e-01, 4.7528820509e-01, 9.2482622270e-02, 7.5419633935e-03])), - np.array([1.2968687054e+00, 9.1901538645e-01, 5.5211228246e-01, 2.5929320267e-01]), - ), - 0.7: ( - mpfutil.normalize( - np.array([4.0619261279e-01, 4.8211835269e-01, 1.0275368140e-01, 8.9353531141e-03])), - np.array([1.3328732728e+00, 9.2538839592e-01, 5.4973718519e-01, 2.5638262685e-01]), - ), - 0.725: ( - mpfutil.normalize( - np.array([3.9117459271e-01, 4.8639603270e-01, 1.1208423379e-01, 1.0345140795e-02])), - np.array([1.3678589145e+00, 9.3142767823e-01, 5.4718845799e-01, 2.5336660763e-01]), - ), - 0.75: ( - mpfutil.normalize( - np.array([3.7883214393e-01, 4.8886181980e-01, 1.2054759480e-01, 1.1758441469e-02])), - np.array([1.4018952939e+00, 9.3712739357e-01, 5.4450683796e-01, 2.5028075100e-01]), - ), - 0.775: ( - mpfutil.normalize( - np.array([3.6858115238e-01, 4.9002935491e-01, 1.2822402555e-01, 1.3165467156e-02])), - np.array([1.4350485124e+00, 9.4249753902e-01, 5.4172689777e-01, 2.4715619082e-01]), - ), - 0.8: ( - mpfutil.normalize( - np.array([3.5999400190e-01, 4.9025638794e-01, 1.3519106673e-01, 1.4558543438e-02])), - np.array([1.4673735837e+00, 9.4755130911e-01, 5.3887419185e-01, 2.4401511404e-01]), - ), - 0.825: ( - mpfutil.normalize( - np.array([3.5275638884e-01, 4.8979354452e-01, 1.4151859016e-01, 1.5931476486e-02])), - np.array([1.4989145047e+00, 9.5229840735e-01, 5.3596493005e-01, 2.4087532891e-01]), - ), - 0.85: ( - mpfutil.normalize( - np.array([3.4661930227e-01, 4.8882885283e-01, 1.4727213167e-01, 1.7279713237e-02])), - np.array([1.5297163559e+00, 9.5675535334e-01, 5.3301353941e-01, 2.3774833319e-01]), - ), - 0.875: ( - mpfutil.normalize( - np.array([3.4139475178e-01, 4.8749519173e-01, 1.5251013392e-01, 1.8599922570e-02])), - np.array([1.5598149092e+00, 9.6093516491e-01, 5.3003111821e-01, 2.3464487823e-01]), - ), - 0.9: ( - mpfutil.normalize( - np.array([3.3693483823e-01, 4.8589217444e-01, 1.5728351471e-01, 1.9889472624e-02])), - np.array([1.5892416384e+00, 9.6484911931e-01, 5.2702460930e-01, 2.3157178352e-01]), - ), - 0.925: ( - mpfutil.normalize( - np.array([3.3311881357e-01, 4.8409461947e-01, 1.6163971512e-01, 2.1146851843e-02])), - np.array([1.6180273453e+00, 9.6851160006e-01, 5.2400218926e-01, 2.2853547559e-01]), - ), - 0.95: ( - mpfutil.normalize( - np.array([3.2985357508e-01, 4.8215861242e-01, 1.6561775424e-01, 2.2370058269e-02])), - np.array([1.6461952712e+00, 9.7193082380e-01, 5.2096478580e-01, 2.2553624276e-01]), - ), - 0.975: ( - mpfutil.normalize( - np.array([3.2705820815e-01, 4.8012679182e-01, 1.6925584733e-01, 2.3559152698e-02])), - np.array([1.6737719728e+00, 9.7512079513e-01, 5.1792043407e-01, 2.2258062250e-01]), - ), - 1.0: ( - mpfutil.normalize( - np.array([3.2466847784e-01, 4.7803210251e-01, 1.7258586760e-01, 2.4713552046e-02])), - np.array([1.7007792632e+00, 9.7809127671e-01, 5.1487137336e-01, 2.1966932363e-01]), - ), - 1.1: ( - mpfutil.normalize( - np.array([3.1821477023e-01, 4.6944264972e-01, 1.8335597803e-01, 2.8986602026e-02])), - np.array([1.8034938165e+00, 9.8796648656e-01, 5.0267797561e-01, 2.0848694531e-01]), - ), - 1.2: ( - mpfutil.normalize( - np.array([3.1519204595e-01, 4.6105964089e-01, 1.9101500145e-01, 3.2733311709e-02])), - np.array([1.8984714894e+00, 9.9502366952e-01, 4.9057127084e-01, 1.9805003387e-01]), - ), - 1.3: ( - mpfutil.normalize( - np.array([3.1432860593e-01, 4.5319073488e-01, 1.9647826914e-01, 3.6002390043e-02])), - np.array([1.9865281111e+00, 9.9969455805e-01, 4.7861874939e-01, 1.8832607467e-01]), - ), - 1.4: ( - mpfutil.normalize( - np.array([3.1489512079e-01, 4.4590727991e-01, 2.0035001581e-01, 3.8847583497e-02])), - np.array([2.0682658389e+00, 1.0023112528e+00, 4.6685531166e-01, 1.7926310694e-01]), - ), - 1.5: ( - mpfutil.normalize( - np.array([3.1644196302e-01, 4.3919486381e-01, 2.0304297828e-01, 4.1320194886e-02])), - np.array([2.1441585351e+00, 1.0031405656e+00, 4.5529969615e-01, 1.7080392320e-01]), - ), - 1.6: ( - mpfutil.normalize( - np.array([3.1868206214e-01, 4.3300475565e-01, 2.0484577176e-01, 4.3467410445e-02])), - np.array([2.2145442254e+00, 1.0023921643e+00, 4.4396221165e-01, 1.6289631422e-01]), - ), - 1.7: ( - mpfutil.normalize( - np.array([3.2141661149e-01, 4.2728356657e-01, 2.0596836872e-01, 4.5331453223e-02])), - np.array([2.2797211146e+00, 1.0002492514e+00, 4.3285136654e-01, 1.5549262413e-01]), - ), - 1.8: ( - mpfutil.normalize( - np.array([3.2450372508e-01, 4.2197913345e-01, 2.0656682602e-01, 4.6950315445e-02])), - np.array([2.3399462224e+00, 9.9687259873e-01, 4.2197602006e-01, 1.4855234090e-01]), - ), - 1.9: ( - mpfutil.normalize( - np.array([3.2782474774e-01, 4.1703843886e-01, 2.0677020639e-01, 4.8366607011e-02])), - np.array([2.3954526149e+00, 9.9243861765e-01, 4.1137384028e-01, 1.4206176868e-01]), - ), - 2.0: ( - mpfutil.normalize( - np.array([3.3128672873e-01, 4.1242201447e-01, 2.0667672958e-01, 4.9614527222e-02])), - np.array([2.4465543239e+00, 9.8711540712e-01, 4.0107011168e-01, 1.3599717559e-01]), - ), - 2.1: ( - mpfutil.normalize( - np.array([3.3480773477e-01, 4.0809304875e-01, 2.0636991637e-01, 5.0729300115e-02])), - np.array([2.4935840432e+00, 9.8108195933e-01, 3.9110438291e-01, 1.3034765947e-01]), - ), - 2.2: ( - mpfutil.normalize( - np.array([3.3824474964e-01, 4.0400630511e-01, 2.0596694378e-01, 5.1782001471e-02])), - np.array([2.5371137448e+00, 9.7470582393e-01, 3.8164575988e-01, 1.2518661787e-01]), - ), - 2.3: ( - mpfutil.normalize( - np.array([3.4154001680e-01, 4.0014083300e-01, 2.0551993522e-01, 5.2799214983e-02])), - np.array([2.5776523647e+00, 9.6817460552e-01, 3.7272302287e-01, 1.2049124896e-01]), - ), - 2.4: ( - mpfutil.normalize( - np.array([3.4462155040e-01, 3.9647171906e-01, 2.0508756596e-01, 5.3819164576e-02])), - np.array([2.6157607079e+00, 9.6172950359e-01, 3.6440804188e-01, 1.1627023797e-01]), - ), - 2.5: ( - mpfutil.normalize( - np.array([3.4743061538e-01, 3.9297757432e-01, 2.0471664490e-01, 5.4875165409e-02])), - np.array([2.6519979386e+00, 9.5559201385e-01, 3.5676044559e-01, 1.1252623389e-01]), - ), - 2.6: ( - mpfutil.normalize( - np.array([3.4967753819e-01, 3.8961227234e-01, 2.0459321850e-01, 5.6116970967e-02])), - np.array([2.6880626669e+00, 9.5066322139e-01, 3.5023361900e-01, 1.0948440882e-01]), - ), - 2.7: ( - mpfutil.normalize( - np.array([3.5177955585e-01, 3.8641279396e-01, 2.0446411884e-01, 5.7343531359e-02])), - np.array([2.7224001800e+00, 9.4582839068e-01, 3.4410642827e-01, 1.0670704905e-01]), - ), - 2.8: ( - mpfutil.normalize( - np.array([3.5319024425e-01, 3.8330557458e-01, 2.0466754604e-01, 5.8836635130e-02])), - np.array([2.7580928560e+00, 9.4277497993e-01, 3.3927833449e-01, 1.0466777695e-01]), - ), - 2.9: ( - mpfutil.normalize( - np.array([3.5378909446e-01, 3.8026742773e-01, 2.0527342180e-01, 6.0670056004e-02])), - np.array([2.7964529303e+00, 9.4202609384e-01, 3.3595005870e-01, 1.0343690690e-01]), - ), - 3.0: ( - mpfutil.normalize( - np.array([3.5433518510e-01, 3.7736713030e-01, 2.0582911847e-01, 6.2468566130e-02])), - np.array([2.8336134738e+00, 9.4133657412e-01, 3.3284256153e-01, 1.0231209515e-01]), - ), - 3.1: ( - mpfutil.normalize( - np.array([3.5408315961e-01, 3.7451250882e-01, 2.0677913224e-01, 6.4625199334e-02])), - np.array([2.8743074645e+00, 9.4315430327e-01, 3.3122820579e-01, 1.0194775480e-01]), - ), - 3.2: ( - mpfutil.normalize( - np.array([3.5394142199e-01, 3.7179082231e-01, 2.0759163388e-01, 6.6676121820e-02])), - np.array([2.9132258233e+00, 9.4461371375e-01, 3.2955559590e-01, 1.0152897182e-01]), - ), - 3.3: ( - mpfutil.normalize( - np.array([3.5305658132e-01, 3.6909677926e-01, 2.0876704539e-01, 6.9079594029e-02])), - np.array([2.9561537781e+00, 9.4863034763e-01, 3.2931721623e-01, 1.0181154872e-01]), - ), - 3.4: ( - mpfutil.normalize( - np.array([3.5192719487e-01, 3.6647515631e-01, 2.1001638061e-01, 7.1581268212e-02])), - np.array([3.0001702025e+00, 9.5363280939e-01, 3.2966495376e-01, 1.0235436427e-01]), - ), - 3.5: ( - mpfutil.normalize( - np.array([3.5060254214e-01, 3.6392322866e-01, 2.1131414032e-01, 7.4160088878e-02])), - np.array([3.0451387321e+00, 9.5950945004e-01, 3.3052314087e-01, 1.0311659760e-01]), - ), - 3.6: ( - mpfutil.normalize( - np.array([3.4912428096e-01, 3.6143851508e-01, 2.1263958405e-01, 7.6797619914e-02])), - np.array([3.0909310508e+00, 9.6616182047e-01, 3.3182712939e-01, 1.0406449740e-01]), - ), - 3.7: ( - mpfutil.normalize( - np.array([3.4752854894e-01, 3.5901873588e-01, 2.1397544966e-01, 7.9477265519e-02])), - np.array([3.1374206755e+00, 9.7350002889e-01, 3.3352041677e-01, 1.0516982075e-01]), - ), - 3.8: ( - mpfutil.normalize( - np.array([3.4584552795e-01, 3.5666183181e-01, 2.1530800221e-01, 8.2184638032e-02])), - np.array([3.1844907916e+00, 9.8144530745e-01, 3.3555479492e-01, 1.0640926144e-01]), - ), - 3.9: ( - mpfutil.normalize( - np.array([3.4362667421e-01, 3.5429634831e-01, 2.1688691349e-01, 8.5190063993e-02])), - np.array([3.2361730736e+00, 9.9188531513e-01, 3.3884215229e-01, 1.0822186544e-01]), - ), - 4.0: ( - mpfutil.normalize( - np.array([3.4184964448e-01, 3.5205834937e-01, 2.1817572509e-01, 8.7916281061e-02])), - np.array([3.2841593573e+00, 1.0008659833e+00, 3.4144548069e-01, 1.0967514839e-01]), - ), - 4.1: ( - mpfutil.normalize( - np.array([3.4005124040e-01, 3.4987797080e-01, 2.1943432458e-01, 9.0636464215e-02])), - np.array([3.3324038190e+00, 1.0102607769e+00, 3.4427950089e-01, 1.1121285477e-01]), - ), - 4.2: ( - mpfutil.normalize( - np.array([3.3780191983e-01, 3.4768053656e-01, 2.2089659664e-01, 9.3620946971e-02])), - np.array([3.3851409074e+00, 1.0220325225e+00, 3.4828215083e-01, 1.1328173586e-01]), - ), - 4.3: ( - mpfutil.normalize( - np.array([3.3601330133e-01, 3.4560948604e-01, 2.2207355647e-01, 9.6303656152e-02])), - np.array([3.4336805116e+00, 1.0321249127e+00, 3.5150240078e-01, 1.1495552057e-01]), - ), - 4.4: ( - mpfutil.normalize( - np.array([3.3424114354e-01, 3.4359102799e-01, 2.2320863101e-01, 9.8959197467e-02])), - np.array([3.4822167800e+00, 1.0424956754e+00, 3.5488215174e-01, 1.1668476113e-01]), - ), - 4.5: ( - mpfutil.normalize( - np.array([3.3249408218e-01, 3.4162359263e-01, 2.2429988001e-01, 1.0158244518e-01])), - np.array([3.5306750929e+00, 1.0531083728e+00, 3.5840316213e-01, 1.1846250846e-01]), - ), - 4.6: ( - mpfutil.normalize( - np.array([3.3077912016e-01, 3.3970523159e-01, 2.2534617890e-01, 1.0416946935e-01])), - np.array([3.5789852971e+00, 1.0639311710e+00, 3.6205013302e-01, 1.2028303071e-01]), - ), - 4.7: ( - mpfutil.normalize( - np.array([3.2910208683e-01, 3.3783437795e-01, 2.2634683555e-01, 1.0671669966e-01])), - np.array([3.6270819206e+00, 1.0749342719e+00, 3.6580889367e-01, 1.2214120589e-01]), - ), - 4.8: ( - mpfutil.normalize( - np.array([3.2746771584e-01, 3.3600943950e-01, 2.2730162001e-01, 1.0922122465e-01])), - np.array([3.6749037219e+00, 1.0860906741e+00, 3.6966693242e-01, 1.2403261813e-01]), - ), - 4.9: ( - mpfutil.normalize( - np.array([3.2624487697e-01, 3.3430423712e-01, 2.2802322150e-01, 1.1142766441e-01])), - np.array([3.7179708914e+00, 1.0953376653e+00, 3.7265528866e-01, 1.2550566759e-01]), - ), - 5.0: ( - mpfutil.normalize( - np.array([3.2469494436e-01, 3.3256655511e-01, 2.2889385343e-01, 1.1384464709e-01])), - np.array([3.7650917059e+00, 1.1067330284e+00, 3.7668264430e-01, 1.2745436857e-01]), - ), - 5.1: ( - mpfutil.normalize( - np.array([3.2354062374e-01, 3.3094533523e-01, 2.2954545188e-01, 1.1596858915e-01])), - np.array([3.8073867328e+00, 1.1161838350e+00, 3.7982690822e-01, 1.2898167857e-01]), - ), - 5.2: ( - mpfutil.normalize( - np.array([3.2208450763e-01, 3.2928880600e-01, 2.3033432453e-01, 1.1829236184e-01])), - np.array([3.8536403099e+00, 1.1277433902e+00, 3.8398921289e-01, 1.3097589146e-01]), - ), - 5.3: ( - mpfutil.normalize( - np.array([3.2100646759e-01, 3.2774571866e-01, 2.3091792559e-01, 1.2032988816e-01])), - np.array([3.8950303683e+00, 1.1373335145e+00, 3.8725931976e-01, 1.3254720884e-01]), - ), - 5.4: ( - mpfutil.normalize( - np.array([3.2027896321e-01, 3.2631406031e-01, 2.3131414720e-01, 1.2209282928e-01])), - np.array([3.9315853626e+00, 1.1449534426e+00, 3.8963891892e-01, 1.3369808741e-01]), - ), - 5.5: ( - mpfutil.normalize( - np.array([3.1926336064e-01, 3.2484343677e-01, 2.3184319316e-01, 1.2405000943e-01])), - np.array([3.9720090223e+00, 1.1546382291e+00, 3.9301530391e-01, 1.3530812924e-01]), - ), - 5.6: ( - mpfutil.normalize( - np.array([3.1858250670e-01, 3.2348118508e-01, 2.3219636882e-01, 1.2573993941e-01])), - np.array([4.0076071518e+00, 1.1623413089e+00, 3.9549599166e-01, 1.3649697272e-01]), - ), - 5.7: ( - mpfutil.normalize( - np.array([3.1792217796e-01, 3.2215146340e-01, 2.3253123511e-01, 1.2739512352e-01])), - np.array([4.0427226616e+00, 1.1700759446e+00, 3.9802186821e-01, 1.3770308584e-01]), - ), - 5.8: ( - mpfutil.normalize( - np.array([3.1728261246e-01, 3.2085314979e-01, 2.3284835496e-01, 1.2901588279e-01])), - np.array([4.0773540484e+00, 1.1778366476e+00, 4.0058956475e-01, 1.3892518114e-01]), - ), - 5.9: ( - mpfutil.normalize( - np.array([3.1666408556e-01, 3.1958517525e-01, 2.3314829004e-01, 1.3060244916e-01])), - np.array([4.1115006498e+00, 1.1856179310e+00, 4.0319577590e-01, 1.4016208250e-01]), - ), - 6.0: ( - mpfutil.normalize( - np.array([3.1633788451e-01, 3.1841919978e-01, 2.3329982319e-01, 1.3194309251e-01])), - np.array([4.1409057150e+00, 1.1914075491e+00, 4.0490227294e-01, 1.4097923007e-01]), - ), - 6.15: ( - mpfutil.normalize( - np.array([3.1599635660e-01, 3.1675771729e-01, 2.3344694983e-01, 1.3379897628e-01])), - np.array([4.1820500944e+00, 1.1991175754e+00, 4.0706101860e-01, 1.4201564715e-01]), - ), - 6.3: ( - mpfutil.normalize( - np.array([3.1566659474e-01, 3.1515443900e-01, 2.3357978807e-01, 1.3559917819e-01])), - np.array([4.2222276786e+00, 1.2068576369e+00, 4.0929437183e-01, 1.4308291759e-01]), - ), - 6.45: ( - mpfutil.normalize( - np.array([3.1559853443e-01, 3.1367747961e-01, 2.3357997276e-01, 1.3714401320e-01])), - np.array([4.2572411244e+00, 1.2126102163e+00, 4.1066204676e-01, 1.4374733334e-01]), - ), - 6.6: ( - mpfutil.normalize( - np.array([3.1577453573e-01, 3.1232227982e-01, 2.3345859196e-01, 1.3844459249e-01])), - np.array([4.2871343136e+00, 1.2163609353e+00, 4.1115794342e-01, 1.4400826897e-01]), - ), - 6.8: ( - mpfutil.normalize( - np.array([3.1615162465e-01, 3.1063558605e-01, 2.3322383402e-01, 1.3998895528e-01])), - np.array([4.3228602263e+00, 1.2200215999e+00, 4.1128396422e-01, 1.4410941934e-01]), - ), - 7.0: ( - mpfutil.normalize( - np.array([3.1673813265e-01, 3.0909602483e-01, 2.3288611110e-01, 1.4127973142e-01])), - np.array([4.3529043052e+00, 1.2216422195e+00, 4.1056166866e-01, 1.4382225816e-01]), - ), - 7.2: ( - mpfutil.normalize( - np.array([3.1775635772e-01, 3.0776545301e-01, 2.3234364820e-01, 1.4213454107e-01])), - np.array([4.3729586302e+00, 1.2191186098e+00, 4.0802622143e-01, 1.4271029990e-01]), - ), - 7.4: ( - mpfutil.normalize( - np.array([3.1896494236e-01, 3.0656756874e-01, 2.3171244184e-01, 1.4275504706e-01])), - np.array([4.3872282722e+00, 1.2144236572e+00, 4.0459360870e-01, 1.4119856069e-01]), - ), - 7.6: ( - mpfutil.normalize( - np.array([3.2036561864e-01, 3.0549662996e-01, 2.3099307346e-01, 1.4314467793e-01])), - np.array([4.3956132281e+00, 1.2074877237e+00, 4.0023879131e-01, 1.3928059243e-01]), - ), - 7.8: ( - mpfutil.normalize( - np.array([3.2196444937e-01, 3.0454758594e-01, 2.3018380480e-01, 1.4330415989e-01])), - np.array([4.3979592256e+00, 1.1982273959e+00, 3.9493334561e-01, 1.3694918085e-01]), - ), - 8.0: ( - mpfutil.normalize( - np.array([3.2401528949e-01, 3.0378290082e-01, 2.2916820523e-01, 1.4303360447e-01])), - np.array([4.3893063163e+00, 1.1843287450e+00, 3.8765576491e-01, 1.3375549788e-01]), - ), - }, - 8: { - 0.5: ( - mpfutil.normalize(np.array([1., 0, 0, 0, 0, 0, 0, 0])), - np.array([1., 9.0907169503e-01, 7.8061774471e-01, 5.6632815891e-01, - 4.0699848261e-01, 2.7428165889e-01, 1.6667433497e-01, 8.1936830104e-02]), - ), - 0.501: ( - mpfutil.normalize(np.array( - [9.8309081062e-01, 1.4087009237e-02, 2.4916406163e-03, 2.4668023661e-04, - 7.3745899515e-05, 7.3233075902e-06, 2.6828887191e-06, 1.0719227828e-07])), - np.array([1.0022362632e+00, 9.1034401393e-01, 7.8116382234e-01, 5.6655450917e-01, - 4.0699848261e-01, 2.7417207769e-01, 1.6654118188e-01, 8.1838662954e-02]), - ), - 0.502: ( - mpfutil.normalize(np.array( - [9.6663028145e-01, 2.7700592068e-02, 4.9990098784e-03, 5.0273560468e-04, - 1.4675979308e-04, 1.5055502427e-05, 5.3505768485e-06, 2.1512466830e-07])), - np.array([1.0044730508e+00, 9.1161557118e-01, 7.8170919194e-01, 5.6678049827e-01, - 4.0699848261e-01, 2.7406275865e-01, 1.6640840035e-01, 8.1740808813e-02]), - ), - 0.503: ( - mpfutil.normalize(np.array( - [9.5060305253e-01, 4.0857769233e-02, 7.5207581369e-03, 7.6773816414e-04, - 2.1921385256e-04, 2.3120652785e-05, 8.0247997539e-06, 3.2263403005e-07])), - np.array([1.0067103642e+00, 9.1288636876e-01, 7.8225385585e-01, 5.6700612752e-01, - 4.0699848261e-01, 2.7395370061e-01, 1.6627598862e-01, 8.1643266061e-02]), - ), - 0.504: ( - mpfutil.normalize(np.array( - [9.3499330605e-01, 5.3577093716e-02, 1.0054040662e-02, 1.0419495572e-03, - 2.9082514962e-04, 3.1712507187e-05, 1.0631785117e-05, 4.4057661952e-07])), - np.array([1.0089482026e+00, 9.1415640864e-01, 7.8279781637e-01, 5.6723139821e-01, - 4.0699848261e-01, 2.7384490244e-01, 1.6614394491e-01, 8.1546033092e-02]), - ), - 0.505: ( - mpfutil.normalize(np.array( - [9.1978663802e-01, 6.5875009211e-02, 1.2596889736e-02, 1.3252955933e-03, - 3.6162779264e-04, 4.0777171354e-05, 1.3184022145e-05, 5.7844986723e-07])), - np.array([1.0111865650e+00, 9.1542569277e-01, 7.8334107580e-01, 5.6745631160e-01, - 4.0699848261e-01, 2.7373636301e-01, 1.6601226749e-01, 8.1449108310e-02]), - ), - 0.506: ( - mpfutil.normalize(np.array( - [9.0497024754e-01, 7.7764895442e-02, 1.5149366619e-02, 1.6168280315e-03, - 4.3206062801e-04, 5.0140682061e-05, 1.5755349674e-05, 7.0570695410e-07])), - np.array([1.0134254504e+00, 9.1669422309e-01, 7.8388363644e-01, 5.6768086897e-01, - 4.0699848261e-01, 2.7362808119e-01, 1.6588095460e-01, 8.1352490130e-02]), - ), - 0.507: ( - mpfutil.normalize(np.array( - [8.9052988922e-01, 8.9263922450e-02, 1.7708231070e-02, 1.9171022956e-03, - 5.0176176043e-04, 5.9970283261e-05, 1.8275338331e-05, 8.4758433427e-07])), - np.array([1.0156648580e+00, 9.1796200154e-01, 7.8442550056e-01, 5.6790507160e-01, - 4.0699848261e-01, 2.7352005585e-01, 1.6575000453e-01, 8.1256176980e-02]), - ), - 0.508: ( - mpfutil.normalize(np.array( - [8.7645376939e-01, 1.0038451549e-01, 2.0273367304e-02, 2.2253623113e-03, - 5.7107404161e-04, 7.0125599132e-05, 2.0793626493e-05, 9.9223547090e-07])), - np.array([1.0179047869e+00, 9.1922903005e-01, 7.8496667043e-01, 5.6812892073e-01, - 4.0699848261e-01, 2.7341228590e-01, 1.6561941557e-01, 8.1160167299e-02]), - ), - 0.509: ( - mpfutil.normalize(np.array( - [8.6272908630e-01, 1.1114182942e-01, 2.2842249206e-02, 2.5419745249e-03, - 6.3968876885e-04, 8.0774130388e-05, 2.3246514064e-05, 1.1511378976e-06])), - np.array([1.0201452361e+00, 9.2049531053e-01, 7.8550714830e-01, 5.6835241762e-01, - 4.0699848261e-01, 2.7330477022e-01, 1.6548918599e-01, 8.1064459535e-02]), - ), - 0.51: ( - mpfutil.normalize(np.array( - [8.4934469629e-01, 1.2154790138e-01, 2.5414362440e-02, 2.8663645915e-03, - 7.0790665304e-04, 9.1766084840e-05, 2.5689755528e-05, 1.3128085204e-06])), - np.array([1.0223862047e+00, 9.2176084490e-01, 7.8604693641e-01, 5.6857556352e-01, - 4.0699848261e-01, 2.7319750772e-01, 1.6535931412e-01, 8.0969052148e-02]), - ), - 0.512: ( - mpfutil.normalize(np.array( - [8.2355249852e-01, 1.4135612329e-01, 3.0563410006e-02, 3.5377502180e-03, - 8.4322152425e-04, 1.1481963468e-04, 3.0522667684e-05, 1.6541427379e-06])), - np.array([1.0268696968e+00, 9.2428968288e-01, 7.8712445227e-01, 5.6902080726e-01, - 4.0699848261e-01, 2.7298373789e-01, 1.6510063678e-01, 8.0779132403e-02]), - ), - 0.514: ( - mpfutil.normalize(np.array( - [7.9899316921e-01, 1.5990380647e-01, 3.5710887171e-02, 4.2384189024e-03, - 9.7711604708e-04, 1.3931148321e-04, 3.5271270163e-05, 2.0194418417e-06])), - np.array([1.0313552560e+00, 9.2681555911e-01, 7.8819923570e-01, 5.6946466178e-01, - 4.0699848261e-01, 2.7277096777e-01, 1.6484337022e-01, 8.0590395959e-02]), - ), - 0.516: ( - mpfutil.normalize(np.array( - [7.7558943163e-01, 1.7727738776e-01, 4.0848854081e-02, 4.9668992165e-03, - 1.1099035731e-03, 1.6516342092e-04, 3.9951716791e-05, 2.4085992607e-06])), - np.array([1.0358428751e+00, 9.2933848856e-01, 7.8927130424e-01, 5.6990713680e-01, - 4.0699848261e-01, 2.7255918881e-01, 1.6458750131e-01, 8.0402830884e-02]), - ), - 0.518: ( - mpfutil.normalize(np.array( - [7.5327013658e-01, 1.9355619523e-01, 4.5970324903e-02, 5.7217595931e-03, - 1.2418748103e-03, 1.9231841549e-04, 4.4569809178e-05, 2.8206549857e-06])), - np.array([1.0403325471e+00, 9.3185848607e-01, 7.9034067523e-01, 5.7034824193e-01, - 4.0699848261e-01, 2.7234839259e-01, 1.6433301709e-01, 8.0216425408e-02]), - ), - 0.52: ( - mpfutil.normalize(np.array( - [7.3196919232e-01, 2.0881407192e-01, 5.1068574515e-02, 6.5018106068e-03, - 1.3732305680e-03, 2.2072873063e-04, 4.9138994140e-05, 3.2523462096e-06])), - np.array([1.0448242651e+00, 9.3437556634e-01, 7.9140736581e-01, 5.7078798666e-01, - 4.0699848261e-01, 2.7213857077e-01, 1.6407990475e-01, 8.0031167928e-02]), - ), - 0.525: ( - mpfutil.normalize(np.array( - [6.8278301465e-01, 2.4292409615e-01, 6.3676783770e-02, 8.5539752337e-03, - 1.7003591408e-03, 2.9698308173e-04, 6.0356445549e-05, 4.4315291824e-06])), - np.array([1.0560624656e+00, 9.4065559947e-01, 7.9406247663e-01, 5.7188145825e-01, - 4.0699848261e-01, 2.7161822553e-01, 1.6345304341e-01, 7.9572971987e-02]), - ), - 0.53: ( - mpfutil.normalize(np.array( - [6.3875507340e-01, 2.7198914397e-01, 7.6033064829e-02, 1.0737137406e-02, - 2.0283922703e-03, 3.8004136368e-04, 7.1417886623e-05, 5.7288694628e-06])), - np.array([1.0673133028e+00, 9.4691771495e-01, 7.9670120149e-01, 5.7296663012e-01, - 4.0699848261e-01, 2.7110379338e-01, 1.6283448496e-01, 7.9121705340e-02]), - ), - 0.535: ( - mpfutil.normalize(np.array( - [5.9918808196e-01, 2.9677542159e-01, 8.8082893011e-02, 1.3034751655e-02, - 2.3599928831e-03, 4.6930724393e-04, 8.2410922273e-05, 7.1407388649e-06])), - np.array([1.0785766716e+00, 9.5316213221e-01, 7.9932379468e-01, 5.7404364267e-01, - 4.0699848261e-01, 2.7059515228e-01, 1.6222404265e-01, 7.8677199293e-02]), - ), - 0.54: ( - mpfutil.normalize(np.array( - [5.6349849628e-01, 3.1791867916e-01, 9.9787377137e-02, 1.5431669887e-02, - 2.6974630300e-03, 5.6423257658e-04, 9.3423425058e-05, 8.6585015208e-06])), - np.array([1.0898524684e+00, 9.5938906599e-01, 8.0193050422e-01, 5.7511263265e-01, - 4.0699848261e-01, 2.7009218384e-01, 1.6162153560e-01, 7.8239290788e-02]), - ), - 0.545: ( - mpfutil.normalize(np.array( - [5.3119503339e-01, 3.3594909797e-01, 1.1111995514e-01, 1.7913961609e-02, - 3.0428527023e-03, 6.6427281862e-04, 1.0455643613e-04, 1.0269938465e-05])), - np.array([1.1011405919e+00, 9.6559872645e-01, 8.0452157205e-01, 5.7617373328e-01, - 4.0699848261e-01, 2.6959477312e-01, 1.6102678861e-01, 7.7807822162e-02]), - ), - 0.55: ( - mpfutil.normalize(np.array( - [5.0186142640e-01, 3.5131169403e-01, 1.2206307201e-01, 2.0469217092e-02, - 3.3977246723e-03, 7.6902605266e-04, 1.1586668614e-04, 1.1973059914e-05])), - np.array([1.1124409422e+00, 9.7179131932e-01, 8.0709723426e-01, 5.7722707438e-01, - 4.0699848261e-01, 2.6910280858e-01, 1.6043963190e-01, 7.7382640919e-02]), - ), - 0.56: ( - mpfutil.normalize(np.array( - [4.5073790392e-01, 3.7547361660e-01, 1.4274718868e-01, 2.5754127142e-02, - 4.1411315564e-03, 9.9111737243e-04, 1.3929326547e-04, 1.5621467014e-05])), - np.array([1.1350779325e+00, 9.8412610389e-01, 8.1220325793e-01, 5.7931098086e-01, - 4.0699848261e-01, 2.6813478777e-01, 1.5928743610e-01, 7.6550555182e-02]), - ), - 0.57: ( - mpfutil.normalize(np.array( - [4.0785785248e-01, 3.9276483697e-01, 1.6182242700e-01, 3.1207637577e-02, - 4.9356424877e-03, 1.2279060332e-03, 1.6413305380e-04, 1.9564395427e-05])), - np.array([1.1577626743e+00, 9.9639498193e-01, 8.1725035376e-01, 5.8136532669e-01, - 4.0699848261e-01, 2.6718729131e-01, 1.5816369065e-01, 7.5741909117e-02]), - ), - 0.58: ( - mpfutil.normalize(np.array( - [3.7153148809e-01, 4.0489676054e-01, 1.7932862668e-01, 3.6765184855e-02, - 5.7860931467e-03, 1.4773544501e-03, 1.9072454716e-04, 2.3767684661e-05])), - np.array([1.1804944278e+00, 1.0085994528e+00, 8.2224021804e-01, 5.8339103895e-01, - 4.0699848261e-01, 2.6625953525e-01, 1.5706721212e-01, 7.4955648308e-02]), - ), - 0.59: ( - mpfutil.normalize(np.array( - [3.4048043140e-01, 4.1312661743e-01, 1.9533862662e-01, 4.2373984655e-02, - 6.6949390946e-03, 1.7378698767e-03, 2.1932542340e-04, 2.8205501710e-05])), - np.array([1.2032724773e+00, 1.0207409564e+00, 8.2717446997e-01, 5.8538900026e-01, - 4.0699848261e-01, 2.6535077843e-01, 1.5599688562e-01, 7.4190782895e-02]), - ), - 0.6: ( - mpfutil.normalize(np.array( - [3.1372329256e-01, 4.1838892696e-01, 2.0994265280e-01, 4.7991121392e-02, - 7.6628074417e-03, 2.0082132442e-03, 2.5012708760e-04, 3.2858515965e-05])), - np.array([1.2260961300e+00, 1.0328208769e+00, 8.3205465644e-01, 5.8736005161e-01, - 4.0699848261e-01, 2.6446031949e-01, 1.5495165976e-01, 7.3446382611e-02]), - ), - 0.625: ( - mpfutil.normalize(np.array( - [2.6099719492e-01, 4.2275873840e-01, 2.4094764459e-01, 6.1858674645e-02, - 1.0333654133e-02, 2.7211565707e-03, 3.3759375082e-04, 4.5342994264e-05])), - np.array([1.2833504226e+00, 1.0627598647e+00, 8.4402813520e-01, 5.9217512270e-01, - 4.0699848261e-01, 2.6230995013e-01, 1.5244203081e-01, 7.1669297996e-02]), - ), - 0.65: ( - mpfutil.normalize(np.array( - [2.2255710120e-01, 4.1964078694e-01, 2.6528960995e-01, 7.5200146771e-02, - 1.3333013103e-02, 3.4795798994e-03, 4.4081020756e-04, 5.8951924994e-05])), - np.array([1.3408756591e+00, 1.0923416226e+00, 8.5569420844e-01, 5.9683849500e-01, - 4.0699848261e-01, 2.6026040245e-01, 1.5006913872e-01, 7.0002434997e-02]), - ), - 0.675: ( - mpfutil.normalize(np.array( - [1.9362331316e-01, 4.1274100861e-01, 2.8430601328e-01, 8.7808324085e-02, - 1.6610115988e-02, 4.2775717437e-03, 5.6000775566e-04, 7.3645380633e-05])), - np.array([1.3986626227e+00, 1.1215839185e+00, 8.6707215232e-01, 6.0136051589e-01, - 4.0699848261e-01, 2.5830333518e-01, 1.4782068946e-01, 6.8435098592e-02]), - ), - 0.7: ( - mpfutil.normalize(np.array( - [1.7126594507e-01, 4.0404691053e-01, 2.9909801838e-01, 9.9583864066e-02, - 2.0109494619e-02, 5.1114689786e-03, 6.9485978320e-04, 8.9438569740e-05])), - np.array([1.4567027409e+00, 1.1505030139e+00, 8.7817937248e-01, 6.0575047484e-01, - 4.0699848261e-01, 2.5643137456e-01, 1.4568589847e-01, 6.6957979020e-02]), - ), - 0.725: ( - mpfutil.normalize(np.array( - [1.5360782153e-01, 3.9464579833e-01, 3.1054479446e-01, 1.1049425719e-01, - 2.3777561692e-02, 5.9786496787e-03, 8.4474182712e-04, 1.0637528435e-04])), - np.array([1.5149880192e+00, 1.1791138408e+00, 8.8903164529e-01, 6.1001674496e-01, - 4.0699848261e-01, 2.5463797213e-01, 1.4365525919e-01, 6.5562930470e-02]), - ), - 0.75: ( - mpfutil.normalize(np.array( - [1.3940075454e-01, 3.8513629592e-01, 3.1933927013e-01, 1.2054767394e-01, - 2.7565695764e-02, 6.8769128425e-03, 1.0088829868e-03, 1.2451388725e-04])), - np.array([1.5735109829e+00, 1.2074301525e+00, 8.9964332072e-01, 6.1416690157e-01, - 4.0699848261e-01, 2.5291728764e-01, 1.4172035363e-01, 6.4242791372e-02]), - ), - 0.775: ( - mpfutil.normalize(np.array( - [1.2778826451e-01, 3.7584414361e-01, 3.2602546393e-01, 1.2977616332e-01, - 3.1431419324e-02, 7.8041692968e-03, 1.1864600849e-03, 1.4391592234e-04])), - np.array([1.6322646278e+00, 1.2354646525e+00, 9.1002749419e-01, 6.1820782198e-01, - 4.0699848261e-01, 2.5126409175e-01, 1.3987369630e-01, 6.2991237363e-02]), - ), - 0.8: ( - mpfutil.normalize(np.array( - [1.1816644753e-01, 3.6693989578e-01, 3.3103051455e-01, 1.3822498822e-01, - 3.5338551735e-02, 8.7583121202e-03, 1.3766491797e-03, 1.6464087983e-04])), - np.array([1.6912423757e+00, 1.2632291058e+00, 9.2019615278e-01, 6.2214576994e-01, - 4.0699848261e-01, 2.4967368486e-01, 1.3810860478e-01, 6.1802660163e-02]), - ), - 0.825: ( - mpfutil.normalize(np.array( - [1.1009934273e-01, 3.5850471394e-01, 3.3469070109e-01, 1.4594581545e-01, - 3.9256854522e-02, 9.7371742451e-03, 1.5786541629e-03, 1.8674385878e-04])), - np.array([1.7504380364e+00, 1.2907344353e+00, 9.3016030047e-01, 6.2598646774e-01, - 4.0699848261e-01, 2.4814182879e-01, 1.3641909160e-01, 6.0672067129e-02]), - ), - 0.85: ( - mpfutil.normalize(np.array( - [1.0326557970e-01, 3.5056829259e-01, 3.3727166915e-01, 1.5299249841e-01, - 4.3161426759e-02, 1.0738540143e-02, 1.7917199410e-03, 2.1027331046e-04])), - np.array([1.8098457733e+00, 1.3179908059e+00, 9.3993006588e-01, 6.2973515774e-01, - 4.0699848261e-01, 2.4666468910e-01, 1.3479977348e-01, 5.9594997451e-02]), - ), - 0.875: ( - mpfutil.normalize(np.array( - [9.7423912383e-02, 3.4313078327e-01, 3.3898426511e-01, 1.5941840593e-01, - 4.7032060966e-02, 1.1760160418e-02, 2.0151406821e-03, 2.3527124063e-04])), - np.array([1.8694600738e+00, 1.3450076975e+00, 9.4951479537e-01, 6.3339665549e-01, - 4.0699848261e-01, 2.4523878609e-01, 1.3324579473e-01, 5.8567451847e-02]), - ), - 0.9: ( - mpfutil.normalize(np.array( - [9.2390237228e-02, 3.3617591771e-01, 3.3999656404e-01, 1.6527489454e-01, - 5.0852565130e-02, 1.2799789712e-02, 2.2482589806e-03, 2.6177265423e-04])), - np.array([1.9292757222e+00, 1.3717939694e+00, 9.5892313392e-01, 6.3697539550e-01, - 4.0699848261e-01, 2.4386095287e-01, 1.3175276218e-01, 5.7585833246e-02]), - ), - 0.925: ( - mpfutil.normalize(np.array( - [8.8021978134e-02, 3.2967878455e-01, 3.4044315871e-01, 1.7061043240e-01, - 5.4610162054e-02, 1.3855212679e-02, 2.4904650251e-03, 2.8980643964e-04])), - np.array([1.9892877762e+00, 1.3983579171e+00, 9.6816309563e-01, 6.4047547107e-01, - 4.0699848261e-01, 2.4252829955e-01, 1.3031668973e-01, 5.6646896507e-02]), - ), - 0.95: ( - mpfutil.normalize(np.array( - [8.4207252704e-02, 3.2361045944e-01, 3.4043231855e-01, 1.7547016950e-01, - 5.8294939748e-02, 1.4924270998e-02, 2.7411934590e-03, 3.1939560181e-04])), - np.array([2.0494915452e+00, 1.4247073224e+00, 9.7724212546e-01, 6.4390066898e-01, - 4.0699848261e-01, 2.4123818220e-01, 1.2893395081e-01, 5.5747705586e-02]), - ), - 0.975: ( - mpfutil.normalize(np.array( - [8.0857203429e-02, 3.1794076520e-01, 3.4005150474e-01, 1.7989578654e-01, - 6.1899379840e-02, 1.6004882352e-02, 2.9999198800e-03, 3.5055801062e-04])), - np.array([2.1098825713e+00, 1.4508494976e+00, 9.8616715342e-01, 6.4725449979e-01, - 4.0699848261e-01, 2.3998817614e-01, 1.2760123754e-01, 5.4885596884e-02]), - ), - 1.0: ( - mpfutil.normalize(np.array( - [7.7900489670e-02, 3.1263985956e-01, 3.3937166369e-01, 1.8392552231e-01, - 6.5417944499e-02, 1.7095056297e-02, 3.2661568475e-03, 3.8330712707e-04])), - np.array([2.1704566125e+00, 1.4767913255e+00, 9.9494464239e-01, 6.5054022453e-01, - 4.0699848261e-01, 2.3877605265e-01, 1.2631552536e-01, 5.4058147753e-02]), - ), - 1.025: ( - mpfutil.normalize(np.array( - [7.6952700357e-02, 3.0727072639e-01, 3.3806252649e-01, 1.8739161187e-01, - 6.8202569142e-02, 1.8172957079e-02, 3.5271339749e-03, 4.1977469774e-04])), - np.array([2.2246007451e+00, 1.4987431217e+00, 1.0027090445e+00, 6.5219634543e-01, - 4.0628168801e-01, 2.3742412651e-01, 1.2510112648e-01, 5.3290608702e-02]), - ), - 1.05: ( - mpfutil.normalize(np.array( - [7.5951490248e-02, 3.0215921072e-01, 3.3682448688e-01, 1.9065840105e-01, - 7.0904528998e-02, 1.9249523626e-02, 3.7943353746e-03, 4.5802310098e-04])), - np.array([2.2794637900e+00, 1.5209999838e+00, 1.0106072058e+00, 6.5389107618e-01, - 4.0558893861e-01, 2.3610532396e-01, 1.2392825164e-01, 5.2558180376e-02]), - ), - 1.075: ( - mpfutil.normalize(np.array( - [7.5135748841e-02, 2.9749043462e-01, 3.3546536698e-01, 1.9358256131e-01, - 7.3462853842e-02, 2.0303293013e-02, 4.0625317696e-03, 4.9720962336e-04])), - np.array([2.3339962749e+00, 1.5427814621e+00, 1.0181712176e+00, 6.5537529000e-01, - 4.0477367330e-01, 2.3472627683e-01, 1.2273485458e-01, 5.1825853206e-02]), - ), - 1.1: ( - mpfutil.normalize(np.array( - [7.4498936511e-02, 2.9320210890e-01, 3.3398176183e-01, 1.9621335234e-01, - 7.5896652559e-02, 2.1336738803e-02, 4.3328997907e-03, 5.3754926888e-04])), - np.array([2.3880979256e+00, 1.5640577695e+00, 1.0254444683e+00, 6.5670779125e-01, - 4.0388057030e-01, 2.3332077650e-01, 1.2154952986e-01, 5.1101902033e-02]), - ), - 1.135: ( - mpfutil.normalize(np.array( - [7.3781546972e-02, 2.8779141761e-01, 3.3188702034e-01, 1.9942751998e-01, - 7.9070960103e-02, 2.2738828766e-02, 4.7078885593e-03, 5.9481767571e-04])), - np.array([2.4635595453e+00, 1.5932581651e+00, 1.0351200522e+00, 6.5822680688e-01, - 4.0243572178e-01, 2.3124840302e-01, 1.1982999923e-01, 5.0086639820e-02]), - ), - 1.165: ( - mpfutil.normalize(np.array( - [7.3255209231e-02, 2.8349066695e-01, 3.3001114663e-01, 2.0193081227e-01, - 8.1698254601e-02, 2.3932622436e-02, 5.0353253230e-03, 6.4596256847e-04])), - np.array([2.5283061848e+00, 1.6180366599e+00, 1.0433454396e+00, 6.5959011063e-01, - 4.0128064060e-01, 2.2954064160e-01, 1.1841431935e-01, 4.9258660768e-02]), - ), - 1.2: ( - mpfutil.normalize(np.array( - [7.3049471911e-02, 2.7922670694e-01, 3.2777409366e-01, 2.0424611111e-01, - 8.4376160956e-02, 2.5219119825e-02, 5.4031544571e-03, 7.0518114381e-04])), - np.array([2.6018171474e+00, 1.6452422737e+00, 1.0517342977e+00, 6.6033371673e-01, - 3.9941275973e-01, 2.2729102882e-01, 1.1665885708e-01, 4.8244979488e-02]), - ), - 1.225: ( - mpfutil.normalize(np.array( - [7.2924434642e-02, 2.7636899184e-01, 3.2617310111e-01, 2.0576519823e-01, - 8.6225063920e-02, 2.6126662214e-02, 5.6678426147e-03, 7.4870541818e-04])), - np.array([2.6544035115e+00, 1.6645338608e+00, 1.0576864274e+00, 6.6090989163e-01, - 3.9813579492e-01, 2.2572975028e-01, 1.1544154676e-01, 4.7552398948e-02]), - ), - 1.25: ( - mpfutil.normalize(np.array( - [7.2878372178e-02, 2.7372507995e-01, 3.2457764447e-01, 2.0711095747e-01, - 8.7973563948e-02, 2.7011369908e-02, 5.9304555740e-03, 7.9255650193e-04])), - np.array([2.7066254051e+00, 1.6834170100e+00, 1.0633867607e+00, 6.6135506613e-01, - 3.9681235914e-01, 2.2415713241e-01, 1.1422435723e-01, 4.6865499137e-02]), - ), - 1.275: ( - mpfutil.normalize(np.array( - [7.2900064724e-02, 2.7126977308e-01, 3.2298627118e-01, 2.0831410943e-01, - 8.9630084705e-02, 2.7871360953e-02, 6.1914750802e-03, 8.3686084817e-04])), - np.array([2.7584803155e+00, 1.7019220657e+00, 1.0688677183e+00, 6.6169363066e-01, - 3.9544448122e-01, 2.2257783114e-01, 1.1301698786e-01, 4.6188041940e-02]), - ), - 1.3: ( - mpfutil.normalize(np.array( - [7.2979341729e-02, 2.6899437999e-01, 3.2142049900e-01, 2.0938091404e-01, - 9.1191523992e-02, 2.8703037565e-02, 6.4489271228e-03, 8.8137655603e-04])), - np.array([2.8099956910e+00, 1.7200600517e+00, 1.0741141966e+00, 6.6190223258e-01, - 3.9401454055e-01, 2.2097674237e-01, 1.1180871339e-01, 4.5518187402e-02]), - ), - 1.333: ( - mpfutil.normalize(np.array( - [7.3183531359e-02, 2.6622654631e-01, 3.1935400320e-01, 2.1060770733e-01, - 9.3134084452e-02, 2.9767976028e-02, 6.7853773974e-03, 9.4077392683e-04])), - np.array([2.8772823086e+00, 1.7433867156e+00, 1.0807222068e+00, 6.6204790231e-01, - 3.9209039282e-01, 2.1886715969e-01, 1.1023231316e-01, 4.4652846949e-02]), - ), - 1.367: ( - mpfutil.normalize(np.array( - [7.3487650222e-02, 2.6368130331e-01, 3.1728975108e-01, 2.1164342964e-01, - 9.4962024290e-02, 3.0809475699e-02, 7.1244411775e-03, 1.0019245866e-03])), - np.array([2.9458863132e+00, 1.7666767337e+00, 1.0870376420e+00, 6.6191991667e-01, - 3.8997933697e-01, 2.1665397927e-01, 1.0860925861e-01, 4.3771069938e-02]), - ), - 1.4: ( - mpfutil.normalize(np.array( - [7.3869856966e-02, 2.6142876206e-01, 3.1531641556e-01, 2.1248342756e-01, - 9.6611627440e-02, 3.1780193876e-02, 7.4479301542e-03, 1.0617863755e-03])), - np.array([3.0116383837e+00, 1.7886160797e+00, 1.0928185226e+00, 6.6164541195e-01, - 3.8788452017e-01, 2.1450781557e-01, 1.0705495422e-01, 4.2939156538e-02]), - ), - 1.433: ( - mpfutil.normalize(np.array( - [7.4328242690e-02, 2.5938643401e-01, 3.1338642759e-01, 2.1316962669e-01, - 9.8136224123e-02, 3.2707434813e-02, 7.7642093822e-03, 1.1214007049e-03])), - np.array([3.0765848120e+00, 1.8098861415e+00, 1.0982259900e+00, 6.6119448681e-01, - 3.8572311817e-01, 2.1234935339e-01, 1.0550848966e-01, 4.2116388403e-02]), - ), - 1.467: ( - mpfutil.normalize(np.array( - [7.4872130853e-02, 2.5746831740e-01, 3.1143663917e-01, 2.1373580324e-01, - 9.9593558843e-02, 3.3625122131e-02, 8.0851260910e-03, 1.1833022694e-03])), - np.array([3.1425944917e+00, 1.8311293882e+00, 1.1034494396e+00, 6.6058831914e-01, - 3.8346504735e-01, 2.1014371612e-01, 1.0394711183e-01, 4.1296657743e-02]), - ), - 1.5: ( - mpfutil.normalize(np.array( - [7.5475995304e-02, 2.5580085499e-01, 3.0959618736e-01, 2.1414861124e-01, - 1.0088311645e-01, 3.4466212022e-02, 8.3863846098e-03, 1.2426380301e-03])), - np.array([3.2056749438e+00, 1.8509983599e+00, 1.1080953316e+00, 6.5978514332e-01, - 3.8117892784e-01, 2.0797174159e-01, 1.0242946309e-01, 4.0505571318e-02]), - ), - 1.55: ( - mpfutil.normalize(np.array( - [7.6502550682e-02, 2.5354880371e-01, 3.0687977586e-01, 2.1456974992e-01, - 1.0265739141e-01, 3.5675680025e-02, 8.8329264259e-03, 1.3331219657e-03])), - np.array([3.2994227931e+00, 1.8798737903e+00, 1.1145220254e+00, 6.5832978192e-01, - 3.7767101636e-01, 2.0472023935e-01, 1.0018979020e-01, 3.9356761547e-02]), - ), - 1.6: ( - mpfutil.normalize(np.array( - [7.7659640625e-02, 2.5160547066e-01, 3.0425539667e-01, 2.1476783169e-01, - 1.0422121322e-01, 3.6802497423e-02, 9.2643975227e-03, 1.4235521771e-03])), - np.array([3.3907987775e+00, 1.9072096944e+00, 1.1201757114e+00, 6.5655601831e-01, - 3.7408162114e-01, 2.0149226610e-01, 9.8004352643e-02, 3.8254670807e-02]), - ), - 1.65: ( - mpfutil.normalize(np.array( - [7.8931342445e-02, 2.4992357670e-01, 3.0171842461e-01, 2.1477754279e-01, - 1.0559902455e-01, 3.7853939452e-02, 9.6820144816e-03, 1.5141349720e-03])), - np.array([3.4797254302e+00, 1.9330494055e+00, 1.1251070506e+00, 6.5450298569e-01, - 3.7043681992e-01, 1.9830296249e-01, 9.5880915235e-02, 3.7202542831e-02]), - ), - 1.7: ( - mpfutil.normalize(np.array( - [8.0316450199e-02, 2.4847646967e-01, 2.9925866462e-01, 2.1461892325e-01, - 1.0680473696e-01, 3.8833724373e-02, 1.0086021036e-02, 1.6050098864e-03])), - np.array([3.5659821303e+00, 1.9573463532e+00, 1.1293146082e+00, 6.5217764373e-01, - 3.6674282182e-01, 1.9515682520e-01, 9.3821522044e-02, 3.6201613818e-02]), - ), - 1.75: ( - mpfutil.normalize(np.array( - [8.1792739839e-02, 2.4721611839e-01, 2.9687267009e-01, 2.1432474982e-01, - 1.0786563839e-01, 3.9752220066e-02, 1.0479168661e-02, 1.6966947385e-03])), - np.array([3.6496319553e+00, 1.9802223376e+00, 1.1328922108e+00, 6.4964332531e-01, - 3.6303886453e-01, 1.9207521841e-01, 9.1836795706e-02, 3.5255752554e-02]), - ), - 1.8: ( - mpfutil.normalize(np.array( - [8.3351702690e-02, 2.4611430142e-01, 2.9455220920e-01, 2.1391460573e-01, - 1.0879847905e-01, 4.0615967685e-02, 1.0863087155e-02, 1.7896470696e-03])), - np.array([3.7306079604e+00, 2.0017116295e+00, 1.1358814976e+00, 6.4693054226e-01, - 3.5934370681e-01, 1.8906845047e-01, 8.9931445695e-02, 3.4366897568e-02]), - ), - 1.8452: ( - mpfutil.normalize(np.array( - [8.4690418058e-02, 2.4499663850e-01, 2.9243076853e-01, 2.1357193298e-01, - 1.0970349732e-01, 4.1464012154e-02, 1.1252890196e-02, 1.8898422607e-03])), - np.array([3.8030355049e+00, 2.0212284065e+00, 1.1390572974e+00, 6.4504670150e-01, - 3.5652992231e-01, 1.8677753901e-01, 8.8522131675e-02, 3.3772539140e-02]), - ), - 1.893: ( - mpfutil.normalize(np.array( - [8.6268380180e-02, 2.4410876298e-01, 2.9029725056e-01, 2.1304087852e-01, - 1.1044271963e-01, 4.2234851844e-02, 1.1620818904e-02, 1.9863373736e-03])), - np.array([3.8760980063e+00, 2.0397751312e+00, 1.1412332930e+00, 6.4239149877e-01, - 3.5318176699e-01, 1.8414199214e-01, 8.6904763730e-02, 3.3059879077e-02]), - ), - 1.942: ( - mpfutil.normalize(np.array( - [8.7934133461e-02, 2.4328016593e-01, 2.8815677809e-01, 2.1243234807e-01, - 1.1112132440e-01, 4.2991916450e-02, 1.1995116441e-02, 2.0882171592e-03])), - np.array([3.9485387128e+00, 2.0576847387e+00, 1.1430695257e+00, 6.3960142547e-01, - 3.4981275648e-01, 1.8153666466e-01, 8.5330623472e-02, 3.2383811008e-02]), - ), - 2.0: ( - mpfutil.normalize(np.array( - [8.9786551909e-02, 2.4210106985e-01, 2.8558295875e-01, 2.1177213071e-01, - 1.1202974355e-01, 4.3990700969e-02, 1.2502640792e-02, 2.2342034665e-03])), - np.array([4.0333353005e+00, 2.0792055401e+00, 1.1460213576e+00, 6.3713731777e-01, - 3.4655738268e-01, 1.7902175966e-01, 8.3868948324e-02, 3.1840955642e-02]), - ), - 2.044: ( - mpfutil.normalize(np.array( - [9.1301957044e-02, 2.4139389885e-01, 2.8371463599e-01, 2.1116596260e-01, - 1.1256844305e-01, 4.4658864885e-02, 1.2856909009e-02, 2.3393285737e-03])), - np.array([4.0944704801e+00, 2.0938548002e+00, 1.1473987381e+00, 6.3482553563e-01, - 3.4386194054e-01, 1.7699616012e-01, 8.2693177817e-02, 3.1382289258e-02]), - ), - 2.1: ( - mpfutil.normalize(np.array( - [9.3065382654e-02, 2.4023047887e-01, 2.8126728442e-01, 2.1048767010e-01, - 1.1340973474e-01, 4.5646379434e-02, 1.3388576606e-02, 2.5044931731e-03])), - np.array([4.1722649038e+00, 2.1134211544e+00, 1.1502793656e+00, 6.3291535695e-01, - 3.4126243863e-01, 1.7503188302e-01, 8.1614864264e-02, 3.1060199486e-02]), - ), - 2.15: ( - mpfutil.normalize(np.array( - [9.4560665893e-02, 2.3907453536e-01, 2.7906268506e-01, 2.0991004396e-01, - 1.1421607553e-01, 4.6592364088e-02, 1.3910229972e-02, 2.6734001343e-03])), - np.array([4.2407794410e+00, 2.1310612983e+00, 1.1533780371e+00, 6.3177072474e-01, - 3.3941875419e-01, 1.7363718122e-01, 8.0895869548e-02, 3.0919426066e-02]), - ), - 2.2: ( - mpfutil.normalize(np.array( - [9.6028518114e-02, 2.3788946552e-01, 2.7686687533e-01, 2.0932325507e-01, - 1.1501950473e-01, 4.7558407790e-02, 1.4457011914e-02, 2.8569615308e-03])), - np.array([4.3078088076e+00, 2.1484065686e+00, 1.1566250205e+00, 6.3089362725e-01, - 3.3783301914e-01, 1.7244418222e-01, 8.0313038134e-02, 3.0855843128e-02]), - ), - 2.25: ( - mpfutil.normalize(np.array( - [9.7470398493e-02, 2.3668081800e-01, 2.7468209559e-01, 2.0872621991e-01, - 1.1581640109e-01, 4.8541160065e-02, 1.5027612045e-02, 3.0552948069e-03])), - np.array([4.3733554241e+00, 2.1654572079e+00, 1.1600060126e+00, 6.3026269866e-01, - 3.3648302866e-01, 1.7143359093e-01, 7.9851408625e-02, 3.0858934533e-02]), - ), - 2.3: ( - mpfutil.normalize(np.array( - [9.8889091627e-02, 2.3545475050e-01, 2.7251038588e-01, 2.0811732384e-01, - 1.1660264907e-01, 4.9536933260e-02, 1.5620444560e-02, 3.2684212673e-03])), - np.array([4.4374176927e+00, 2.1821996431e+00, 1.1635002124e+00, 6.2985372602e-01, - 3.3534616120e-01, 1.7058665943e-01, 7.9497067166e-02, 3.0919392113e-02]), - ), - 2.35: ( - mpfutil.normalize(np.array( - [1.0028533292e-01, 2.3421673847e-01, 2.7035450425e-01, 2.0749554501e-01, - 1.1737529596e-01, 5.0542325903e-02, 1.6233944343e-02, 3.4963131519e-03])), - np.array([4.5000020129e+00, 2.1986350310e+00, 1.1670911245e+00, 6.2964634637e-01, - 3.3440248537e-01, 1.6988674041e-01, 7.9237722794e-02, 3.1029137360e-02]), - ), - 2.4: ( - mpfutil.normalize(np.array( - [1.0166002192e-01, 2.3296987374e-01, 2.6821637086e-01, 2.0686191605e-01, - 1.1813193981e-01, 5.1554435824e-02, 1.6866573154e-02, 3.7388686379e-03])), - np.array([4.5611299918e+00, 2.2147646425e+00, 1.1707691398e+00, 6.2962335682e-01, - 3.3363435510e-01, 1.6931896869e-01, 7.9062329480e-02, 3.1180991276e-02]), - ), - 2.45: ( - mpfutil.normalize(np.array( - [1.0301380953e-01, 2.3171819539e-01, 2.6609819914e-01, 2.0621611457e-01, - 1.1887039516e-01, 5.2570453009e-02, 1.7516850242e-02, 3.9959829505e-03])), - np.array([4.6208228082e+00, 2.2305898027e+00, 1.1745208042e+00, 6.2976766696e-01, - 3.3302562238e-01, 1.6887020365e-01, 7.8961503011e-02, 3.1369070707e-02]), - ), - 2.5: ( - mpfutil.normalize(np.array( - [1.0434667717e-01, 2.3046410911e-01, 2.6400142175e-01, 2.0555926083e-01, - 1.1958933776e-01, 5.3588202632e-02, 1.8183451128e-02, 4.2675396129e-03])), - np.array([4.6791134755e+00, 2.2461191985e+00, 1.1783397851e+00, 6.3006658874e-01, - 3.3256311146e-01, 1.6852944693e-01, 7.8927226712e-02, 3.1588428026e-02]), - ), - 2.55: ( - mpfutil.normalize(np.array( - [1.0566067184e-01, 2.2921303327e-01, 2.6192774405e-01, 2.0488960950e-01, - 1.2028636846e-01, 5.4604654370e-02, 1.8864649102e-02, 4.5532694101e-03])), - np.array([4.7360096495e+00, 2.2613349155e+00, 1.1822032180e+00, 6.3050051791e-01, - 3.3223081812e-01, 1.6828481685e-01, 7.8951447628e-02, 3.1834441327e-02]), - ), - 2.6: ( - mpfutil.normalize(np.array( - [1.0663061559e-01, 2.2745849245e-01, 2.5964517834e-01, 2.0436860610e-01, - 1.2130437859e-01, 5.5923644186e-02, 1.9740214167e-02, 4.9288705730e-03])), - np.array([4.7970061271e+00, 2.2803015813e+00, 1.1890070616e+00, 6.3307311259e-01, - 3.3338463451e-01, 1.6903101454e-01, 7.9605696314e-02, 3.2456364347e-02]), - ), - 2.65: ( - mpfutil.normalize(np.array( - [1.0790557056e-01, 2.2621886621e-01, 2.5761869750e-01, 2.0367233833e-01, - 1.2195285249e-01, 5.6934476185e-02, 2.0451267627e-02, 5.2459310922e-03])), - np.array([4.8513257908e+00, 2.2949999186e+00, 1.1929968082e+00, 6.3377421628e-01, - 3.3329580065e-01, 1.6895976256e-01, 7.9732103854e-02, 3.2746172362e-02]), - ), - 2.7: ( - mpfutil.normalize(np.array( - [1.0916065205e-01, 2.2498629638e-01, 2.5561907907e-01, 2.0296754679e-01, - 1.2257791626e-01, 5.7938884730e-02, 2.1173112678e-02, 5.5765120366e-03])), - np.array([4.9043712073e+00, 2.3094181753e+00, 1.1970194368e+00, 6.3458200588e-01, - 3.3330793029e-01, 1.6896082545e-01, 7.9900586569e-02, 3.3053266578e-02]), - ), - 2.75: ( - mpfutil.normalize(np.array( - [1.1008444908e-01, 2.2327735208e-01, 2.5340785630e-01, 2.0238676034e-01, - 1.2350278648e-01, 5.9239679667e-02, 2.2095060009e-02, 6.0060560484e-03])), - np.array([4.9617054917e+00, 2.3276741994e+00, 1.2040126981e+00, 6.3752656632e-01, - 3.3478996844e-01, 1.6993382919e-01, 8.0683460395e-02, 3.3725853600e-02]), - ), - 2.8: ( - mpfutil.normalize(np.array( - [1.1130256368e-01, 2.2206871038e-01, 2.5146278144e-01, 2.0166019969e-01, - 1.2407490850e-01, 6.0226911651e-02, 2.2837963326e-02, 6.3659613192e-03])), - np.array([5.0123574283e+00, 2.3415859718e+00, 1.2081076172e+00, 6.3853444592e-01, - 3.3498546135e-01, 1.7006268720e-01, 8.0923919786e-02, 3.4060491447e-02]), - ), - 2.85: ( - mpfutil.normalize(np.array( - [1.1220074338e-01, 2.2040401673e-01, 2.4930663492e-01, 2.0104169873e-01, - 1.2493053945e-01, 6.1503513038e-02, 2.3782704733e-02, 6.8301490227e-03])), - np.array([5.0673631713e+00, 2.3593519738e+00, 1.2151688819e+00, 6.4166361061e-01, - 3.3663257234e-01, 1.7114745588e-01, 8.1767184064e-02, 3.4753608041e-02]), - ), - 2.9: ( - mpfutil.normalize(np.array( - [1.1308692527e-01, 2.1876701579e-01, 2.4717939847e-01, 2.0040280720e-01, - 1.2574853032e-01, 6.2765221618e-02, 2.4737699842e-02, 7.3124014987e-03])), - np.array([5.1212891084e+00, 2.3768980296e+00, 1.2222642154e+00, 6.4487729962e-01, - 3.3835302626e-01, 1.7227977187e-01, 8.2634402516e-02, 3.5453447574e-02]), - ), - 2.95: ( - mpfutil.normalize(np.array( - [1.1425537366e-01, 2.1760738089e-01, 2.4532216450e-01, 1.9964738926e-01, - 1.2623603774e-01, 6.3714296381e-02, 2.5502860637e-02, 7.7144969257e-03])), - np.array([5.1685573726e+00, 2.3900620555e+00, 1.2264140867e+00, 6.4611837614e-01, - 3.3876464178e-01, 1.7255731719e-01, 8.2956029128e-02, 3.5816189477e-02]), - ), - 3.0: ( - mpfutil.normalize(np.array( - [1.1512004551e-01, 2.1602259133e-01, 2.4325646168e-01, 1.9897828593e-01, - 1.2698319194e-01, 6.4940948770e-02, 2.6470606715e-02, 8.2278681347e-03])), - np.array([5.2202744700e+00, 2.4070910785e+00, 1.2335136248e+00, 6.4945365553e-01, - 3.4060010334e-01, 1.7376728150e-01, 8.3863632609e-02, 3.6527471879e-02]), - ), - 3.05: ( - mpfutil.normalize(np.array( - [1.1597603466e-01, 2.1446903244e-01, 2.4122283549e-01, 1.9829262197e-01, - 1.2769190000e-01, 6.6147103833e-02, 2.7443267828e-02, 8.7572037825e-03])), - np.array([5.2709410709e+00, 2.4238712126e+00, 1.2406093168e+00, 6.5283996166e-01, - 3.4248386203e-01, 1.7500877486e-01, 8.4786459362e-02, 3.7241868618e-02]), - ), - 3.1: ( - mpfutil.normalize(np.array( - [1.1682327549e-01, 2.1294678841e-01, 2.3922179756e-01, 1.9759210514e-01, - 1.2836298994e-01, 6.7331851687e-02, 2.8419406682e-02, 9.3017850934e-03])), - np.array([5.3205873987e+00, 2.4404050386e+00, 1.2476961525e+00, 6.5627146080e-01, - 3.4441115790e-01, 1.7627853183e-01, 8.5722658846e-02, 3.7958672257e-02]), - ), - 3.15: ( - mpfutil.normalize(np.array( - [1.1791931900e-01, 2.1185627919e-01, 2.3748530895e-01, 1.9681559060e-01, - 1.2874711731e-01, 6.8221432892e-02, 2.9200460543e-02, 9.7544915254e-03])), - np.array([5.3639723435e+00, 2.4527440049e+00, 1.2519177111e+00, 6.5777262574e-01, - 3.4505590302e-01, 1.7670998382e-01, 8.6123277006e-02, 3.8345208220e-02]), - ), - 3.2: ( - mpfutil.normalize(np.array( - [1.1874159838e-01, 2.1038521507e-01, 2.3554812627e-01, 1.9609760389e-01, - 1.2935521717e-01, 6.9366597752e-02, 3.0180228075e-02, 1.0325413402e-02])), - np.array([5.4117238630e+00, 2.4688349151e+00, 1.2589983675e+00, 6.6129517461e-01, - 3.4706894200e-01, 1.7803544069e-01, 8.7086292086e-02, 3.9067758155e-02]), - ), - 3.25: ( - mpfutil.normalize(np.array( - [1.1955490715e-01, 2.0894489038e-01, 2.3364407234e-01, 1.9536925413e-01, - 1.2992860305e-01, 7.0488734974e-02, 3.1160007761e-02, 1.0909530209e-02])), - np.array([5.4585475724e+00, 2.4846949676e+00, 1.2660619546e+00, 6.6485043614e-01, - 3.4911482362e-01, 1.7938190215e-01, 8.8058614860e-02, 3.9791255852e-02]), - ), - 3.3: ( - mpfutil.normalize(np.array( - [1.2035894604e-01, 2.0753504215e-01, 2.3177347256e-01, 1.9463180441e-01, - 1.3046827304e-01, 7.1587533143e-02, 3.2138792937e-02, 1.1506135726e-02])), - np.array([5.5044785584e+00, 2.5003304861e+00, 1.2731057625e+00, 6.6843459236e-01, - 3.5119081978e-01, 1.8074764735e-01, 8.9039347408e-02, 4.0515387900e-02]), - ), - 3.35: ( - mpfutil.normalize(np.array( - [1.2115419218e-01, 2.0615509329e-01, 2.2993537788e-01, 1.9388697112e-01, - 1.3097564604e-01, 7.2662851165e-02, 3.3115439886e-02, 1.2114428439e-02])), - np.array([5.5495324864e+00, 2.5157417771e+00, 1.2801283448e+00, 6.7204544968e-01, - 3.5329429773e-01, 1.8213062261e-01, 9.0027338718e-02, 4.1239870850e-02]), - ), - 3.4: ( - mpfutil.normalize(np.array( - [1.2194010012e-01, 2.0480463414e-01, 2.2813066047e-01, 1.9313573237e-01, - 1.3145163960e-01, 7.3714337390e-02, 3.4089146264e-02, 1.2733749648e-02])), - np.array([5.5937486897e+00, 2.5309389611e+00, 1.2871275756e+00, 6.7567916387e-01, - 3.5542268340e-01, 1.8352951659e-01, 9.1022085410e-02, 4.1964551740e-02]), - ), - 3.45: ( - mpfutil.normalize(np.array( - [1.2271679126e-01, 2.0348313844e-01, 2.2635874198e-01, 1.9237936436e-01, - 1.3189749153e-01, 7.4742082356e-02, 3.5059028879e-02, 1.3363361193e-02])), - np.array([5.6371509832e+00, 2.5459257457e+00, 1.2941024292e+00, 6.7933376534e-01, - 3.5757415065e-01, 1.8494301007e-01, 9.2022843779e-02, 4.2689251287e-02]), - ), - 3.5: ( - mpfutil.normalize(np.array( - [1.2348421641e-01, 2.0219007573e-01, 2.2461937267e-01, 1.9161899885e-01, - 1.3231440362e-01, 7.5746083321e-02, 3.6024285201e-02, 1.4002564198e-02])), - np.array([5.6797665163e+00, 2.5607077696e+00, 1.3010519273e+00, 6.8300709141e-01, - 3.5974680430e-01, 1.8636987221e-01, 9.3029040552e-02, 4.3413847119e-02]), - ), - 3.55: ( - mpfutil.normalize(np.array( - [1.2424269973e-01, 2.0092510736e-01, 2.2291222875e-01, 1.9085554443e-01, - 1.3270338934e-01, 7.6726270318e-02, 3.6984124800e-02, 1.4650635282e-02])), - np.array([5.7216133126e+00, 2.5752862911e+00, 1.3079731089e+00, 6.8669613262e-01, - 3.6193835992e-01, 1.8780871178e-01, 9.4040057156e-02, 4.4138218449e-02]), - ), - 3.6: ( - mpfutil.normalize(np.array( - [1.2499190637e-01, 1.9968730238e-01, 2.2123681947e-01, 1.9009007645e-01, - 1.3306576562e-01, 7.7683187886e-02, 3.7938031883e-02, 1.5306909928e-02])), - np.array([5.7627256419e+00, 2.5896713573e+00, 1.3148684086e+00, 6.9040093047e-01, - 3.6414857894e-01, 1.8925924496e-01, 9.5055564287e-02, 4.4862260965e-02]), - ), - 3.65: ( - mpfutil.normalize(np.array( - [1.2573176524e-01, 1.9847597296e-01, 2.1959281696e-01, 1.8932379628e-01, - 1.3340263296e-01, 7.8616914330e-02, 3.8885387505e-02, 1.5970713767e-02])), - np.array([5.8031293506e+00, 2.6038701081e+00, 1.3217383735e+00, 6.9411996181e-01, - 3.6637578974e-01, 1.9072050284e-01, 9.6075158409e-02, 4.5585915607e-02]), - ), - 3.7: ( - mpfutil.normalize(np.array( - [1.2646305704e-01, 1.9729132514e-01, 2.1797991077e-01, 1.8855679385e-01, - 1.3371477378e-01, 7.9527322669e-02, 3.9825450713e-02, 1.6641366042e-02])), - np.array([5.8428298493e+00, 2.6178766687e+00, 1.3285757967e+00, 6.9784880404e-01, - 3.6861737786e-01, 1.9219092415e-01, 9.7098285314e-02, 4.6309119118e-02]), - ), - 3.75: ( - mpfutil.normalize(np.array( - [1.2718511680e-01, 1.9613214067e-01, 2.1639764460e-01, 1.8779049881e-01, - 1.3400350178e-01, 8.0415031735e-02, 4.0757848240e-02, 1.7318217370e-02])), - np.array([5.8818675998e+00, 2.6317064173e+00, 1.3353863421e+00, 7.0158859854e-01, - 3.7087332294e-01, 1.9367040029e-01, 9.8124694603e-02, 4.7031789460e-02]), - ), - 3.8: ( - mpfutil.normalize(np.array( - [1.2789854715e-01, 1.9499807335e-01, 2.1484533224e-01, 1.8702524657e-01, - 1.3426977985e-01, 8.1280206817e-02, 4.1682119953e-02, 1.8000694078e-02])), - np.array([5.9202509209e+00, 2.6453574729e+00, 1.3421669494e+00, 7.0533735334e-01, - 3.7314230377e-01, 1.9515818943e-01, 9.9154193319e-02, 4.7753948042e-02]), - ), - 3.85: ( - mpfutil.normalize(np.array( - [1.2842780769e-01, 1.9361534521e-01, 2.1313202385e-01, 1.8625537170e-01, - 1.3465832295e-01, 8.2323607041e-02, 4.2776120770e-02, 1.8811400791e-02])), - np.array([5.9625696784e+00, 2.6622666168e+00, 1.3514490204e+00, 7.1085520547e-01, - 3.7660561559e-01, 1.9742654944e-01, 1.0067859995e-01, 4.8777439146e-02]), - ), - 3.9: ( - mpfutil.normalize(np.array( - [1.2912868611e-01, 1.9253705987e-01, 2.1164152228e-01, 1.8549125228e-01, - 1.3487685924e-01, 8.3140071516e-02, 4.3680864584e-02, 1.9503684108e-02])), - np.array([5.9996657443e+00, 2.6755408084e+00, 1.3581472204e+00, 7.1460423716e-01, - 3.7888871564e-01, 1.9892341102e-01, 1.0171000379e-01, 4.9496722902e-02]), - ), - 3.95: ( - mpfutil.normalize(np.array( - [1.2982086138e-01, 1.9148203685e-01, 2.1017996197e-01, 1.8473051034e-01, - 1.3507592954e-01, 8.3934873426e-02, 4.4576152096e-02, 2.0199674414e-02])), - np.array([6.0361736941e+00, 2.6886533526e+00, 1.3648152211e+00, 7.1835811171e-01, - 3.8118123021e-01, 2.0042645010e-01, 1.0274359203e-01, 5.0215392865e-02]), - ), - 4.0: ( - mpfutil.normalize(np.array( - [1.3050468325e-01, 1.9044951688e-01, 2.0874610747e-01, 1.8397333270e-01, - 1.3525663822e-01, 8.4708868885e-02, 4.5461937689e-02, 2.0898914905e-02])), - np.array([6.0721045194e+00, 2.7016061284e+00, 1.3714545129e+00, 7.2211812610e-01, - 3.8348404027e-01, 2.0193608252e-01, 1.0377947788e-01, 5.0933501664e-02]), - ), - 4.1: ( - mpfutil.normalize(np.array( - [1.3184825454e-01, 1.8845150671e-01, 2.0596125095e-01, 1.8247090737e-01, - 1.3556538880e-01, 8.6194805579e-02, 4.7203140322e-02, 2.2304745724e-02])), - np.array([6.1422983042e+00, 2.7270382985e+00, 1.3846329619e+00, 7.2964380225e-01, - 3.8811250115e-01, 2.0497049539e-01, 1.0585627643e-01, 5.2367819802e-02]), - ), - 4.2: ( - mpfutil.normalize(np.array( - [1.3301513494e-01, 1.8631141709e-01, 2.0311396247e-01, 1.8096505533e-01, - 1.3591730858e-01, 8.7769491409e-02, 4.9066063914e-02, 2.3841566259e-02])), - np.array([6.2146435327e+00, 2.7550582491e+00, 1.4000653887e+00, 7.3884173113e-01, - 3.9388979167e-01, 2.0875869812e-01, 1.0840933854e-01, 5.4090081334e-02]), - ), - 4.3: ( - mpfutil.normalize(np.array( - [1.3430431293e-01, 1.8449156691e-01, 2.0054270431e-01, 1.7949957768e-01, - 1.3609498671e-01, 8.9091499317e-02, 5.0718536359e-02, 2.5256815789e-02])), - np.array([6.2806592599e+00, 2.7792508586e+00, 1.4129488972e+00, 7.4635052679e-01, - 3.9855381588e-01, 2.1181821820e-01, 1.1049171900e-01, 5.5516568425e-02]), - ), - 4.4: ( - mpfutil.normalize(np.array( - [1.3556410874e-01, 1.8274805186e-01, 1.9806987002e-01, 1.7805825360e-01, - 1.3622032776e-01, 9.0340786925e-02, 5.2327199805e-02, 2.6671401292e-02])), - np.array([6.3447965811e+00, 2.8029000648e+00, 1.4257072603e+00, 7.5385502234e-01, - 4.0323626364e-01, 2.1489091680e-01, 1.1257796183e-01, 5.6940376121e-02]), - ), - 4.5: ( - mpfutil.normalize(np.array( - [1.3679594823e-01, 1.8107788964e-01, 1.9569083703e-01, 1.7664236084e-01, - 1.3629837769e-01, 9.1520760271e-02, 5.3891466965e-02, 2.8082359335e-02])), - np.array([6.4071473021e+00, 2.8260258462e+00, 1.4383381698e+00, 7.6135081278e-01, - 4.0793349707e-01, 2.1797470518e-01, 1.1466727484e-01, 5.8361643837e-02]), - ), - 4.6: ( - mpfutil.normalize(np.array( - [1.3800108456e-01, 1.7947759431e-01, 1.9340174659e-01, 1.7525331410e-01, - 1.3633363569e-01, 9.2634597524e-02, 5.5411072816e-02, 2.9486954411e-02])), - np.array([6.4678088704e+00, 2.8486506844e+00, 1.4508423168e+00, 7.6883372946e-01, - 4.1264180393e-01, 2.2106747670e-01, 1.1675891531e-01, 5.9780434911e-02]), - ), - 4.7: ( - mpfutil.normalize(np.array( - [1.3918127169e-01, 1.7794501050e-01, 1.9119871681e-01, 1.7389082026e-01, - 1.3633023670e-01, 9.3685459716e-02, 5.6885876583e-02, 3.0882607748e-02])), - np.array([6.5268491379e+00, 2.8707816551e+00, 1.4632109827e+00, 7.7629656263e-01, - 4.1735689916e-01, 2.2416663456e-01, 1.1885159191e-01, 6.1196311899e-02]), - ), - 4.8: ( - mpfutil.normalize(np.array( - [1.4033716723e-01, 1.7647598941e-01, 1.8907813850e-01, 1.7255669992e-01, - 1.3629206863e-01, 9.4676513678e-02, 5.8316176201e-02, 3.2267246434e-02])), - np.array([6.5843648296e+00, 2.8924513040e+00, 1.4754531109e+00, 7.8373918855e-01, - 4.2207679887e-01, 2.2727115737e-01, 1.2094518804e-01, 6.2609681569e-02]), - ), - 4.9: ( - mpfutil.normalize(np.array( - [1.4157125392e-01, 1.7522838121e-01, 1.8716966594e-01, 1.7128851232e-01, - 1.3616253195e-01, 9.5491108829e-02, 5.9568398417e-02, 3.3520147420e-02])), - np.array([6.6366410301e+00, 2.9108913690e+00, 1.4854449240e+00, 7.8964888920e-01, - 4.2577282028e-01, 2.2970121174e-01, 1.2260089345e-01, 6.3746109656e-02]), - ), - 5.0: ( - mpfutil.normalize(np.array( - [1.4267733172e-01, 1.7387113205e-01, 1.8519870416e-01, 1.7001196979e-01, - 1.3606987408e-01, 9.6377393452e-02, 6.0915101900e-02, 3.4878492853e-02])), - np.array([6.6913690235e+00, 2.9317311731e+00, 1.4974620436e+00, 7.9706074560e-01, - 4.3050678750e-01, 2.3281907675e-01, 1.2469829023e-01, 6.5156094457e-02]), - ), - 5.1: ( - mpfutil.normalize(np.array( - [1.4385476816e-01, 1.7271653749e-01, 1.8342561300e-01, 1.6880344270e-01, - 1.3590032299e-01, 9.7102754325e-02, 6.2092023240e-02, 3.6104538095e-02])), - np.array([6.7410850083e+00, 2.9494624758e+00, 1.5072860987e+00, 8.0296697109e-01, - 4.3423169553e-01, 2.3527145513e-01, 1.2636356164e-01, 6.6292616115e-02]), - ), - 5.2: ( - mpfutil.normalize(np.array( - [1.4491621392e-01, 1.7146216671e-01, 1.8159203583e-01, 1.6758451266e-01, - 1.3576315676e-01, 9.7893935269e-02, 6.3357091479e-02, 3.7430887362e-02])), - np.array([6.7932718323e+00, 2.9695276563e+00, 1.5190746582e+00, 8.1033303262e-01, - 4.3896750675e-01, 2.3839580208e-01, 1.2846199275e-01, 6.7698849263e-02]), - ), - 5.3: ( - mpfutil.normalize(np.array( - [1.4604271678e-01, 1.7039389811e-01, 1.7994226969e-01, 1.6643347488e-01, - 1.3556178198e-01, 9.8539227951e-02, 6.4461234748e-02, 3.8625395859e-02])), - np.array([6.8406552922e+00, 2.9865933131e+00, 1.5287243865e+00, 8.1622345086e-01, - 4.4271157529e-01, 2.4086452609e-01, 1.3013382272e-01, 6.8834656271e-02]), - ), - 5.4: ( - mpfutil.normalize(np.array( - [1.4722412010e-01, 1.6949366122e-01, 1.7846318982e-01, 1.6535139442e-01, - 1.3530810191e-01, 9.9053342293e-02, 6.5414944160e-02, 3.9691246071e-02])), - np.array([6.8834009876e+00, 3.0007780552e+00, 1.5363056270e+00, 8.2067796810e-01, - 4.4548686862e-01, 2.4269140525e-01, 1.3138700961e-01, 6.9704202530e-02]), - ), - 5.5: ( - mpfutil.normalize(np.array( - [1.4829637111e-01, 1.6849368091e-01, 1.7691939685e-01, 1.6425685308e-01, - 1.3508817294e-01, 9.9635017023e-02, 6.6456106191e-02, 4.0854401893e-02])), - np.array([6.9287658032e+00, 3.0173188729e+00, 1.5458415305e+00, 8.2657886747e-01, - 4.4926293579e-01, 2.4518413066e-01, 1.3307059556e-01, 7.0842731983e-02]), - ), - 5.6: ( - mpfutil.normalize(np.array( - [1.4942165771e-01, 1.6765121824e-01, 1.7553685105e-01, 1.6322992755e-01, - 1.3482308095e-01, 1.0009495176e-01, 6.7353118274e-02, 4.1889194464e-02])), - np.array([6.9695675064e+00, 3.0310137040e+00, 1.5533208942e+00, 8.3104642017e-01, - 4.5207105571e-01, 2.4703541710e-01, 1.3433575994e-01, 7.1715139327e-02]), - ), - 5.7: ( - mpfutil.normalize(np.array( - [1.5051799722e-01, 1.6683328855e-01, 1.7419762841e-01, 1.6222969356e-01, - 1.3455525798e-01, 1.0053213477e-01, 6.8224693805e-02, 4.2909305717e-02])), - np.array([7.0095276785e+00, 3.0445210399e+00, 1.5607750677e+00, 8.3553303435e-01, - 4.5490191327e-01, 2.4890256681e-01, 1.3560934858e-01, 7.2590558062e-02]), - ), - 5.8: ( - mpfutil.normalize(np.array( - [1.5158792725e-01, 1.6604005473e-01, 1.7289989309e-01, 1.6125472516e-01, - 1.3428484264e-01, 1.0094748236e-01, 6.9071232434e-02, 4.3913842346e-02])), - np.array([7.0486425661e+00, 3.0578222685e+00, 1.5681891740e+00, 8.4002957177e-01, - 4.5775019026e-01, 2.5078226758e-01, 1.3688916818e-01, 7.3467835723e-02]), - ), - 5.9: ( - mpfutil.normalize(np.array( - [1.5263279582e-01, 1.6527061144e-01, 1.7164248857e-01, 1.6030482148e-01, - 1.3401191311e-01, 1.0134133028e-01, 6.9893311577e-02, 4.4902727718e-02])), - np.array([7.0869326327e+00, 3.0709231376e+00, 1.5755609077e+00, 8.4453112005e-01, - 4.6061222379e-01, 2.5267309951e-01, 1.3817487885e-01, 7.4346947017e-02]), - ), - 6.0: ( - mpfutil.normalize(np.array( - [1.5371713898e-01, 1.6463100011e-01, 1.7052243811e-01, 1.5941948048e-01, - 1.3371154242e-01, 1.0163821022e-01, 7.0589635597e-02, 4.5770554090e-02])), - np.array([7.1210215057e+00, 3.0813908883e+00, 1.5809925533e+00, 8.4766436584e-01, - 4.6254287480e-01, 2.5394296737e-01, 1.3905432302e-01, 7.4966879418e-02]), - ), - 6.15: ( - mpfutil.normalize(np.array( - [1.5532199260e-01, 1.6375229101e-01, 1.6895141678e-01, 1.5815510695e-01, - 1.3325313744e-01, 1.0202071497e-01, 7.1550315388e-02, 4.6995024867e-02])), - np.array([7.1691377015e+00, 3.0956215831e+00, 1.5881735110e+00, 8.5171833318e-01, - 4.6500955889e-01, 2.5556294912e-01, 1.4018414564e-01, 7.5772718074e-02]), - ), - 6.3: ( - mpfutil.normalize(np.array( - [1.5686603179e-01, 1.6290521063e-01, 1.6744659939e-01, 1.5693971192e-01, - 1.3280243138e-01, 1.0237690308e-01, 7.2472727947e-02, 4.8190383863e-02])), - np.array([7.2157687319e+00, 3.1095524059e+00, 1.5953244210e+00, 8.5581352855e-01, - 4.6752173591e-01, 2.5721456866e-01, 1.4133137743e-01, 7.6585577304e-02]), - ), - 6.45: ( - mpfutil.normalize(np.array( - [1.5840765803e-01, 1.6218346059e-01, 1.6609566205e-01, 1.5581108055e-01, - 1.3233942017e-01, 1.0264089699e-01, 7.3265856087e-02, 4.9255965537e-02])), - np.array([7.2576107819e+00, 3.1208221651e+00, 1.6005851953e+00, 8.5859482922e-01, - 4.6914535235e-01, 2.5827458822e-01, 1.4208894217e-01, 7.7146971445e-02]), - ), - 6.6: ( - mpfutil.normalize(np.array( - [1.5994306187e-01, 1.6157753990e-01, 1.6488859503e-01, 1.5476619314e-01, - 1.3186912057e-01, 1.0282192383e-01, 7.3937826800e-02, 5.0195738872e-02])), - np.array([7.2947660041e+00, 3.1294755712e+00, 1.6039725302e+00, 8.6006890360e-01, - 4.6988388751e-01, 2.5874522320e-01, 1.4245864655e-01, 7.7458236566e-02]), - ), - 6.75: ( - mpfutil.normalize(np.array( - [1.6141910294e-01, 1.6098767282e-01, 1.6372898648e-01, 1.5376188229e-01, - 1.3141305425e-01, 1.0299136793e-01, 7.4585487837e-02, 5.1112445437e-02])), - np.array([7.3306902820e+00, 3.1379271793e+00, 1.6073614821e+00, 8.6159121219e-01, - 4.7066869743e-01, 2.5924786377e-01, 1.4284641961e-01, 7.7777411618e-02]), - ), - 7.0: ( - mpfutil.normalize(np.array( - [1.6386662083e-01, 1.6024541094e-01, 1.6209433699e-01, 1.5226568372e-01, - 1.3064631376e-01, 1.0310845885e-01, 7.5410830021e-02, 5.2362344881e-02])), - np.array([7.3800711240e+00, 3.1460298465e+00, 1.6086521566e+00, 8.6106440904e-01, - 4.6989867480e-01, 2.5870391810e-01, 1.4258610237e-01, 7.7726701690e-02]), - ), - 7.25: ( - mpfutil.normalize(np.array( - [1.6625045892e-01, 1.5970370924e-01, 1.6073156290e-01, 1.5094575689e-01, - 1.2989424076e-01, 1.0309376513e-01, 7.6017051878e-02, 5.3363454282e-02])), - np.array([7.4195180989e+00, 3.1488186642e+00, 1.6061580697e+00, 8.5790733902e-01, - 4.6735872259e-01, 2.5698876004e-01, 1.4155871604e-01, 7.7182431792e-02]), - ), - 7.5: ( - mpfutil.normalize(np.array( - [1.6857173662e-01, 1.5935270491e-01, 1.5962467309e-01, 1.4979276842e-01, - 1.2916003978e-01, 1.0295932451e-01, 7.6416704059e-02, 5.4122048623e-02])), - np.array([7.4490163730e+00, 3.1461993534e+00, 1.5997758459e+00, 8.5204749480e-01, - 4.6300731542e-01, 2.5408172029e-01, 1.3975555905e-01, 7.6141918540e-02]), - ), - 7.75: ( - mpfutil.normalize(np.array( - [1.7083435424e-01, 1.5918958798e-01, 1.5876424273e-01, 1.4879934469e-01, - 1.2844440326e-01, 1.0271170289e-01, 7.6616636570e-02, 5.4639727635e-02])), - np.array([7.4684071981e+00, 3.1379717367e+00, 1.5893304661e+00, 8.4337165764e-01, - 4.5678037973e-01, 2.4994977872e-01, 1.3716165853e-01, 7.4599734900e-02]), - ), - 8.0: ( - mpfutil.normalize(np.array( - [1.7304825187e-01, 1.5921892478e-01, 1.5814669452e-01, 1.4796027259e-01, - 1.2774431184e-01, 1.0235044003e-01, 7.6617131731e-02, 5.4913972642e-02])), - np.array([7.4772419651e+00, 3.1237557818e+00, 1.5745497550e+00, 8.3171078484e-01, - 4.4858426833e-01, 2.4454696379e-01, 1.3375692440e-01, 7.2548667757e-02]), - ), - } - } - } - - @classmethod - def _checkengine(cls, engine): - if engine not in mpfobj.Model.ENGINES: - raise ValueError("Unknown {:s} rendering engine {:s}".format(type(cls), engine)) - - def is_gaussian(self): - return self.profile == "sersic" and self.parameters["n_ser"].value == 0.5 - - def is_gaussian_mixture(self): - return True - - def get_parameters(self, free=True, fixed=True): - return super().get_parameters(free=free, fixed=fixed) + \ - [value for value in self.parameters.values() if - (value.fixed and fixed) or (not value.fixed and free)] - - def get_profiles(self, flux_by_band, engine, cen_x, cen_y, params=None, engineopts=None, - get_derivatives=True): - self._checkengine(engine) - if params is None: - params = {} - - for band in flux_by_band.keys(): - if band not in self.fluxes_dict: - raise ValueError( - f"Asked for EllipticalComponent (profile={self.profile}, name={self.name})" - f"model for band={band} not in bands with fluxes {self.fluxes_dict}") - - profile = {param.name: param.value for param in - self.parameters.values()} - is_sersic = self.profile == "sersic" - slope = profile["n_ser"] if is_sersic else None - if is_sersic: - param_n_ser = self.parameters["n_ser"] - if slope is None: - raise RuntimeError("Can't get multigaussian profiles for profile {}".format(profile)) - can_return_gauss = self.profile == "sersic" and (slope <= 0.5) and ( - not self.do_return_all_profiles or self.parameters["n_ser"].fixed) and not get_derivatives - if can_return_gauss: - weights = [1] - reffs = [1] - profiles = [{}] - else: - slope_log10 = np.log10(slope) - if slope in MultiGaussianApproximationComponent.weights[self.profile][self.order]: - weights, reffs = MultiGaussianApproximationComponent.weights[self.profile][self.order][slope] - negatives = weights < 0 - else: - weights = np.array([f[0](slope_log10) for f in self.weight_splines]) - reffs = np.array([f[0](slope_log10) for f in self.reff_splines]) - negatives = weights < 0 - weights[negatives] = 0 - if not all([(0 <= w <= 1) and (s > 0) for w, s in zip(weights, reffs)]): - raise RuntimeError('Weights {} not all >= 0 and <= 1 and/or reffs {} not all >=0 for slope ' - '{:.4e} log10(slope)={:.4e}'.format(weights, reffs, slope, slope_log10)) - weights, total = mpfutil.normalize(weights, return_sum=True) - if get_derivatives: - try: - # Return dy/dn_ser, not dy/d(log10(n_ser)) - dconst = 1.0/(slope*ln10) - dweights = np.array([f[1](slope_log10)*dconst for f in self.weight_splines]) - dreffs = np.array([f[1](slope_log10)*dconst for f in self.reff_splines]) - except Exception as e: - raise e - dweights[negatives] = 0 - profiles = [{} for _ in range(self.order)] - - profile_base = super().get_profiles(flux_by_band, engine, cen_x, cen_y, params, engineopts)[0] - skip_covar = engineopts is not None and engineopts.get("get_profile_skip_covar", False) - reff, axrat, ang = (np.Inf, None, None) if skip_covar else self.params_ellipse.make_major() - profile_base['n_ser'] = 0.5 - profile_base['can_do_fit_leastsq'] = True - for band in flux_by_band.keys(): - flux = self.fluxes_dict[band].value - profile = profile_base.copy() - if isinstance(self.fluxes_dict[band], g2f.ProperFractionParameterD): - fluxratio = np.float64(flux) - if not 0 <= fluxratio <= 1: - raise ValueError("flux ratio not 0 <= {} <= 1".format(fluxratio)) - flux *= flux_by_band[band] - flux_by_band[band] *= (1.0-fluxratio) - if not skip_covar: - profile["axrat"] = axrat - profile["ang"] = ang - if not 0 < profile["axrat"] <= 1: - if profile["axrat"] <= __class__.axrat_min: - profile["axrat"] = __class__.axrat_min - else: - raise ValueError("axrat {} ! >0 and <=1".format(profile["axrat"])) - - cens = {"cen_x": cen_x, "cen_y": cen_y} - for key, value in cens.items(): - if key in profile: - profile[key] += value - else: - profile[key] = np.float64(value) - if engine == "galsim": - axrat = profile["axrat"] - axrat_sqrt = np.sqrt(axrat) - gsparams = mpfobj.get_gsparams(engineopts) - elif engine == "libprofit": - profile["profile"] = "sersic" - else: - raise ValueError("Unimplemented rendering engine {:s}".format(engine)) - profile["pointsource"] = False - profile["resolved"] = True - profile["param_flux"] = self.fluxes_dict[band] - if is_sersic: - profile["param_n_ser"] = param_n_ser - - for subcomp, (weight, size_reff) in enumerate(zip(weights, reffs)): - size = size_reff*gauss2d.M_HWHM_SIGMA - profile_sub = copy.copy(profile) - # Not needed right now. Maybe later? - #profile_sub["fluxfrac_sub"] = weight - profile_sub["sigma_factor"] = size - profile_sub["sigma_x"] *= size - profile_sub["sigma_y"] *= size - fluxsub = weight*flux - if not fluxsub >= 0: - print(np.array([f[0](slope_log10) for f in self.weight_splines])) - print(weights) - print(reffs) - print(weight, reff, slope_log10, profile_sub) - raise RuntimeError('MGA {} fluxsub={} !>=0'.format(self.profile, fluxsub)) - if engine == "galsim": - profile_gs = gs.Gaussian(flux=weight*flux, half_light_radius=reff*size*axrat_sqrt, - gsparams=gsparams) - profile_sub.update({ - "profile_gs": profile_gs, - "shear": gs.Shear(q=axrat, beta=profile["ang"]*gs.degrees), - "offset": gs.PositionD(profile["cen_x"], profile["cen_y"]), - }) - elif engine == "libprofit": - profile_sub["flux"] = weight * flux - profile_sub["re"] = reff * size - if get_derivatives: - profile_sub["dweight_dn"] = dweights[subcomp] - profile_sub["dreff_dn"] = dreffs[subcomp] - profiles[subcomp][band] = profile_sub - - return profiles - - @classmethod - def _checkparameters(cls, parameters, profile, order): - if profile != "sersic" and order not in MultiGaussianApproximationComponent.weights[profile]: - raise ValueError("{} profile={} order={} not in supported {}".format( - cls.__name__, profile, order, MultiGaussianApproximationComponent.weights[profile].keys())) - - mandatory = {param: False for param in mpfobj.EllipticalParametricComponent.mandatory[profile]} - name_params_needed = mandatory.keys() - name_params = [param.name for param in parameters] - errors = [] - # Not as efficient as early exit if true but oh well - if len(name_params) > len(set(name_params)): - errors.append("Parameters array not unique") - # Check if these parameters are known (in mandatory) - for param in parameters: - if isinstance(param, g2f.IntegralParameterD): - errors.append("Param {:s} is {:s}, not {:s}".format(param.name, type(g2f.IntegralParameterD), - type(mpfobj.Parameter))) - if param.name in name_params_needed: - mandatory[param.name] = True - if param.name == "n_ser": - n_ser = param.value - n_sers = [x for x in MultiGaussianApproximationComponent.weights[profile][order]] - n_ser_min = min(n_sers) - n_ser_max = max(n_sers) - if n_ser < n_ser_min or n_ser > n_ser_max: - raise RuntimeError("Asked for Multigaussian_sersic with n={} not {}= 0 for value in weight_values]) - weight_valuestouse = weight_values[is_weight] - for j, (splines, ext, order_spline) in enumerate( - [(self.weight_splines, 'zeros', 3), (self.reff_splines, 'const', 5)]): - spline = spinterp.InterpolatedUnivariateSpline( - indices[is_weight], [values[j][i] for values in weight_valuestouse], - ext=ext, k=order_spline) - splines.append((spline, spline.derivative())) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..0a6f516 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,106 @@ +[build-system] +requires = ["setuptools", "lsst-versions >= 1.3.0"] +build-backend = "setuptools.build_meta" + +[project] +name = "lsst-multiprofit" +description = "Astronomical image and source model fitting code." +license = {file = "LICENSE"} +readme = "README.rst" +authors = [ + {name="Rubin Observatory Data Management", email="dm-admin@lists.lsst.org"}, +] +classifiers = [ + "Intended Audience :: Science/Research", + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Topic :: Scientific/Engineering :: Astronomy", +] +keywords = [ + "astronomy", + "astrophysics", + "fitting", + "lsst", + "models", + "modeling", +] +requires-python = ">=3.10.0" +dependencies = [ + "astropy", + "gauss2d", + "gauss2dfit", + "lsst-pex-config", + "lsst-utils", + "importlib_resources", + "matplotlib", + "numpy", + "pydantic", + "scipy", +] +dynamic = ["version"] + +[project.urls] +"Homepage" = "https://github.com/lsst-dm/multiprofit" + +[project.optional-dependencies] +galsim = ["galsim"] +test = [ + "pytest", +] + +[tool.setuptools.packages.find] +where = ["python"] + +[tool.setuptools.dynamic] +version = { attr = "lsst_versions.get_lsst_version" } + +[tool.black] +line-length = 110 +target-version = ["py311"] + +[tool.isort] +profile = "black" +line_length = 110 + +[tool.ruff] +exclude = [ + "__init__.py", + "examples/fithsc.py", + "examples/test_utils.py", + "tests/*.py", +] +ignore = [ + "N802", + "N803", + "N806", + "N812", + "N815", + "N816", + "N999", + "D107", + "D105", + "D102", + "D104", + "D100", + "D200", + "D205", + "D400", +] +line-length = 110 +select = [ + "E", # pycodestyle + "F", # pycodestyle + "N", # pep8-naming + "W", # pycodestyle + "D", # pydocstyle +] +target-version = "py311" + +[tool.ruff.pycodestyle] +max-doc-length = 79 + +[tool.ruff.pydocstyle] +convention = "numpy" diff --git a/python/lsst/__init__.py b/python/lsst/__init__.py new file mode 100644 index 0000000..d38efca --- /dev/null +++ b/python/lsst/__init__.py @@ -0,0 +1,4 @@ +import pkgutil + +__path__ = pkgutil.extend_path(__path__, __name__) + diff --git a/python/lsst/multiprofit/__init__.py b/python/lsst/multiprofit/__init__.py new file mode 100644 index 0000000..d38efca --- /dev/null +++ b/python/lsst/multiprofit/__init__.py @@ -0,0 +1,4 @@ +import pkgutil + +__path__ = pkgutil.extend_path(__path__, __name__) + diff --git a/multiprofit/asinhstretchsigned.py b/python/lsst/multiprofit/asinhstretchsigned.py similarity index 85% rename from multiprofit/asinhstretchsigned.py rename to python/lsst/multiprofit/asinhstretchsigned.py index 341fec8..47158f7 100644 --- a/multiprofit/asinhstretchsigned.py +++ b/python/lsst/multiprofit/asinhstretchsigned.py @@ -31,9 +31,8 @@ def _prepare(values, clip=True, out=None): Prepare the data by optionally clipping and copying, and return the array that should be subsequently used for in-place calculations. """ - if clip: - return np.clip(values, 0., 1., out=out) + return np.clip(values, 0.0, 1.0, out=out) else: if out is None: return np.array(values, copy=True) @@ -59,7 +58,7 @@ class AsinhStretchSigned(apvis.BaseStretch): to logarithmic behavior, expressed as a fraction of the normalized image. Must be in the range between 0 and 1. Default is 0.1 - """ + """ # noqa: W505 def __init__(self, a=0.1): super().__init__() @@ -74,14 +73,14 @@ def __call__(self, values, clip=True, out=None): np.abs(values, out=values) np.true_divide(values, self.a, out=values) np.arcsinh(values, out=values) - np.true_divide(values, np.arcsinh(1. / self.a), out=values) - np.true_divide(1. + signs*values, 2., out=values) + np.true_divide(values, np.arcsinh(1.0 / self.a), out=values) + np.true_divide(1.0 + signs * values, 2.0, out=values) return values @property def inverse(self): """A stretch object that performs the inverse operation.""" - return SinhStretchSigned(a=1. / np.arcsinh(1. / self.a)) + return SinhStretchSigned(a=1.0 / np.arcsinh(1.0 / self.a)) class SinhStretchSigned(apvis.BaseStretch): @@ -99,24 +98,24 @@ class SinhStretchSigned(apvis.BaseStretch): The ``a`` parameter used in the above formula. Default is 1/3. """ - def __init__(self, a=1. / 3.): + def __init__(self, a=1.0 / 3.0): super().__init__() self.a = a -# [docs] + # [docs] def __call__(self, values, clip=True, out=None): values = _prepare(values, clip=clip, out=out) - values *= 2. - values -= 1. + values *= 2.0 + values -= 1.0 np.true_divide(values, self.a, out=values) np.sinh(values, out=values) - np.true_divide(values, np.sinh(1. / self.a), out=values) - values += 1. - values /= 2. + np.true_divide(values, np.sinh(1.0 / self.a), out=values) + values += 1.0 + values /= 2.0 return values @property def inverse(self): """A stretch object that performs the inverse operation.""" - return AsinhStretchSigned(a=1. / np.sinh(1. / self.a)) + return AsinhStretchSigned(a=1.0 / np.sinh(1.0 / self.a)) diff --git a/multiprofit/componentconfig.py b/python/lsst/multiprofit/componentconfig.py similarity index 69% rename from multiprofit/componentconfig.py rename to python/lsst/multiprofit/componentconfig.py index 1703043..09ea913 100644 --- a/multiprofit/componentconfig.py +++ b/python/lsst/multiprofit/componentconfig.py @@ -27,17 +27,31 @@ from .transforms import transforms_ref parameter_names = { - g2f.CentroidXParameterD: 'cen_x', - g2f.CentroidYParameterD: 'cen_y', - g2f.ReffXParameterD: 'reff_x', - g2f.ReffYParameterD: 'reff_y', - g2f.RhoParameterD: 'rho', - g2f.SigmaXParameterD: 'sigma_x', - g2f.SigmaYParameterD: 'sigma_y', + g2f.CentroidXParameterD: "cen_x", + g2f.CentroidYParameterD: "cen_y", + g2f.ReffXParameterD: "reff_x", + g2f.ReffYParameterD: "reff_y", + g2f.RhoParameterD: "rho", + g2f.SigmaXParameterD: "sigma_x", + g2f.SigmaYParameterD: "sigma_y", } def init_component(component: g2f.Component, **kwargs): + """Initialize a component with parameter name-value pairs. + + Parameters + ---------- + component + The component to initialize. + kwargs + Additional keyword arguments. + + Notes + ----- + kwargs keywords should be a value in parameter_names and values should be + valid for initializing the parameter of that type. + """ for parameter in set(component.parameters()): if kwarg := parameter_names.get(parameter.__class__): if value := kwargs.get(kwarg): @@ -46,18 +60,14 @@ def init_component(component: g2f.Component, **kwargs): class ParameterConfig(pexConfig.Config): """Basic configuration for all parameters.""" + fixed = pexConfig.Field[bool](default=False, doc="Whether parameter is fixed or not (free)") value_initial = pexConfig.Field[float](default=0, doc="Initial value") -class SersicIndexConfig(ParameterConfig): - """Specific configuration for a Sersic index parameter.""" - def setDefaults(self): - self.value_initial = 0.5 - - class EllipticalComponentConfig(ShapePriorConfig): """Config for an elliptically-symmetric component""" + rho = pexConfig.ConfigField[ParameterConfig](doc="Rho parameter config") size = pexConfig.ConfigField[ParameterConfig](doc="Size parameter config") @@ -94,32 +104,37 @@ def make_component( class GaussianConfig(EllipticalComponentConfig): """Configuration for a gauss2d.fit Gaussian component""" + def make_component( self, centroid: g2f.CentroidParameters, channels: list[g2f.Channel], ) -> g2f.Component: - transform_flux = transforms_ref['log10'] - transform_size = transforms_ref['log10'] - transform_rho = transforms_ref['logit_rho'] + transform_flux = transforms_ref["log10"] + transform_size = transforms_ref["log10"] + transform_rho = transforms_ref["logit_rho"] ellipse = g2f.GaussianParametricEllipse( - sigma_x=g2f.SigmaXParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), - sigma_y=g2f.SigmaYParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), + sigma_x=g2f.SigmaXParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), + sigma_y=g2f.SigmaYParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), rho=g2f.RhoParameterD(self.rho.value_initial, transform=transform_rho, fixed=self.rho.fixed), ) prior = self.get_shape_prior(ellipse) return g2f.GaussianComponent( centroid=centroid, ellipse=ellipse, - integral=g2f.LinearIntegralModel([ - (channel, g2f.IntegralParameterD(1.0, transform=transform_flux)) - for channel in channels - ]), + integral=g2f.LinearIntegralModel( + [(channel, g2f.IntegralParameterD(1.0, transform=transform_flux)) for channel in channels] + ), ), ([] if prior is None else [prior]) class SersicIndexConfig(ParameterConfig): """Specific configuration for a Sersic index parameter.""" + def setDefaults(self): self.value_initial = 0.5 @@ -132,29 +147,37 @@ class SersicConfig(EllipticalComponentConfig): make_component will return a `gauss2d.fit.GaussianComponent` if the Sersic index is fixed at 0.5, or a `gauss2d.fit.SersicMixComponent` otherwise. """ + sersicindex = pexConfig.ConfigField[SersicIndexConfig](doc="Sersic index config") def make_component( - self, - centroid: g2f.CentroidParameters, - channels: list[g2f.Channel], - label_integral: str | None = None + self, centroid: g2f.CentroidParameters, channels: list[g2f.Channel], label_integral: str | None = None ) -> tuple[g2f.Component, list[g2f.Prior]]: is_gaussian = self.sersicindex.value_initial == 0.5 and self.sersicindex.fixed if label_integral is None: label_integral = f"{'Gaussian' if is_gaussian else 'Sersic'} {{channel.name}}-band" - transform_flux = transforms_ref['log10'] - transform_size = transforms_ref['log10'] - transform_rho = transforms_ref['logit_rho'] + transform_flux = transforms_ref["log10"] + transform_size = transforms_ref["log10"] + transform_rho = transforms_ref["logit_rho"] integral = g2f.LinearIntegralModel( - [(channel, g2f.IntegralParameterD(1.0, transform=transform_flux, - label=label_integral.format(channel=channel))) - for channel in channels] + [ + ( + channel, + g2f.IntegralParameterD( + 1.0, transform=transform_flux, label=label_integral.format(channel=channel) + ), + ) + for channel in channels + ] ) if is_gaussian: ellipse = g2f.GaussianParametricEllipse( - sigma_x=g2f.SigmaXParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), - sigma_y=g2f.SigmaYParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), + sigma_x=g2f.SigmaXParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), + sigma_y=g2f.SigmaYParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), rho=g2f.RhoParameterD(self.rho.value_initial, transform=transform_rho, fixed=self.rho.fixed), ) component = g2f.GaussianComponent( @@ -164,8 +187,12 @@ def make_component( ) else: ellipse = g2f.SersicParametricEllipse( - size_x=g2f.ReffXParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), - size_y=g2f.ReffYParameterD(self.size.value_initial, transform=transform_size, fixed=self.size.fixed), + size_x=g2f.ReffXParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), + size_y=g2f.ReffYParameterD( + self.size.value_initial, transform=transform_size, fixed=self.size.fixed + ), rho=g2f.RhoParameterD(self.rho.value_initial, transform=transform_rho, fixed=self.rho.fixed), ) component = g2f.SersicMixComponent( @@ -175,7 +202,7 @@ def make_component( sersicindex=g2f.SersicMixComponentIndexParameterD( value=self.sersicindex.value_initial, fixed=self.sersicindex.fixed, - transform=transforms_ref['logit_sersic'] if not self.sersicindex.fixed else None, + transform=transforms_ref["logit_sersic"] if not self.sersicindex.fixed else None, ), ) prior = self.get_shape_prior(ellipse) diff --git a/multiprofit/config.py b/python/lsst/multiprofit/config.py similarity index 88% rename from multiprofit/config.py rename to python/lsst/multiprofit/config.py index d4cafa5..19d8f70 100644 --- a/multiprofit/config.py +++ b/python/lsst/multiprofit/config.py @@ -27,6 +27,15 @@ def set_config_from_dict( config: pexConfig.Config | pexConfig.dictField.Dict | pexConfig.configDictField.ConfigDict | dict, overrides: dict[str, Any], ): + """Set `lsst.pex.config` params from a dict. + + Parameters + ---------- + config + A config, dictField or configDictField object. + overrides + A dict of key-value pairs to override in the config. + """ is_config_dict = hasattr(config, "__getitem__") if is_config_dict: keys = tuple(config.keys()) diff --git a/multiprofit/fit_bootstrap_model.py b/python/lsst/multiprofit/fit_bootstrap_model.py similarity index 69% rename from multiprofit/fit_bootstrap_model.py rename to python/lsst/multiprofit/fit_bootstrap_model.py index 2209381..c26671d 100644 --- a/multiprofit/fit_bootstrap_model.py +++ b/python/lsst/multiprofit/fit_bootstrap_model.py @@ -24,56 +24,66 @@ from functools import cached_property import gauss2d as g2 import gauss2d.fit as g2f -from multiprofit.config import set_config_from_dict -from multiprofit.componentconfig import init_component -from multiprofit.fit_psf import CatalogExposurePsfABC, CatalogPsfFitterConfig -from multiprofit.fit_source import ( - CatalogExposureSourcesABC, CatalogSourceFitterABC, CatalogSourceFitterConfig, +import numpy + +from .config import set_config_from_dict +from .componentconfig import init_component +from .fit_psf import CatalogExposurePsfABC, CatalogPsfFitterConfig +from .fit_source import ( + CatalogExposureSourcesABC, + CatalogSourceFitterABC, + CatalogSourceFitterConfig, ) -from multiprofit.utils import get_params_uniq +from .utils import get_params_uniq import numpy as np from typing import Any, Mapping @dataclass(kw_only=True, frozen=True) class SourceCatalogBootstrap: + """Config class for a bootstrap source catalog fitter.""" + n_cols_img: int = 25 n_rows_img: int = 25 n_sources: int = 1 @cached_property def catalog(self): - catalog = astropy.table.Table({'id': np.arange(self.n_sources)}) + catalog = astropy.table.Table({"id": np.arange(self.n_sources)}) return catalog @dataclass(kw_only=True, frozen=True) class CatalogExposurePsfBootstrap(CatalogExposurePsfABC, SourceCatalogBootstrap): + """Dataclass for a PSF-convolved bootstrap fitter.""" + noise: float = 1e-4 sigma_x: float sigma_y: float rho: float @cached_property - def centroid(self): - centroid = g2.Centroid(x=self.n_cols_img/2, y=self.n_rows_img/2) + def centroid(self) -> g2.Centroid: + centroid = g2.Centroid(x=self.n_cols_img / 2, y=self.n_rows_img / 2) return centroid @cached_property - def ellipse(self): - g2.Centroid(x=self.n_cols_img/2, y=self.n_rows_img/2) + def ellipse(self) -> g2.Ellipse: + g2.Centroid(x=self.n_cols_img / 2, y=self.n_rows_img / 2) ellipse = g2.Ellipse(g2.EllipseValues(self.sigma_x, self.sigma_y, self.rho)) return ellipse @cached_property - def image(self): + def image(self) -> numpy.ndarray: image = g2.make_gaussians_pixel_D( - g2.ConvolvedGaussians([ - g2.ConvolvedGaussian( - g2.Gaussian(centroid=self.centroid, ellipse=self.ellipse), - g2.Gaussian(), - ) - ]), + g2.ConvolvedGaussians( + [ + g2.ConvolvedGaussian( + g2.Gaussian(centroid=self.centroid, ellipse=self.ellipse), + g2.Gaussian(), + ) + ] + ), n_rows=self.n_rows_img, n_cols=self.n_cols_img, ).data @@ -82,13 +92,15 @@ def image(self): def get_catalog(self) -> astropy.table.Table: return self.catalog - def get_psf_image(self, source): - rng = np.random.default_rng(source['id']) - return self.image + 1e-4*rng.standard_normal(self.image.shape) + def get_psf_image(self, source) -> numpy.ndarray: + rng = np.random.default_rng(source["id"]) + return self.image + 1e-4 * rng.standard_normal(self.image.shape) @dataclass(kw_only=True, frozen=True) class CatalogExposureSourcesBootstrap(CatalogExposureSourcesABC, SourceCatalogBootstrap): + """A CatalogExposure for bootstrap fitting of source catalogs.""" + channel: g2f.Channel = g2f.Channel.NONE config_fit: CatalogSourceFitterConfig model_source: g2f.Source @@ -99,7 +111,7 @@ def get_catalog(self) -> astropy.table.Table: return self.catalog def get_psfmodel(self, params: Mapping[str, Any]) -> g2f.PsfModel: - return self.config_fit_psf.rebuild_psfmodel(self.table_psf_fits[params['id']]) + return self.config_fit_psf.rebuild_psfmodel(self.table_psf_fits[params["id"]]) def get_source_observation(self, source: Mapping[str, Any]) -> g2f.Observation: n_rows = self.n_rows_img + self.n_buffer_img @@ -113,10 +125,10 @@ def get_source_observation(self, source: Mapping[str, Any]) -> g2f.Observation: return obs def __post_init__(self): - config_dict = self.table_psf_fits.meta['config'] + config_dict = self.table_psf_fits.meta["config"] config = CatalogPsfFitterConfig() set_config_from_dict(config, config_dict) - object.__setattr__(self, 'config_fit_psf', config) + object.__setattr__(self, "config_fit_psf", config) @dataclass(kw_only=True) @@ -127,6 +139,7 @@ class CatalogSourceFitterBootstrap(CatalogSourceFitterABC): each row. The resulting catalog can be used to examine performance and statistics of the best-fit parameters. """ + background: float = 1e2 flux: float = 1e4 sigma_x: float @@ -136,36 +149,39 @@ class CatalogSourceFitterBootstrap(CatalogSourceFitterABC): def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float) -> tuple[float, float]: return float(cen_x), float(cen_y) - def initialize_model(self, model: g2f.Model, source: g2f.Source, - limits_x: g2f.LimitsD = None, limits_y: g2f.LimitsD = None): + def initialize_model( + self, model: g2f.Model, source: g2f.Source, limits_x: g2f.LimitsD = None, limits_y: g2f.LimitsD = None + ): comp1, comp2 = model.sources[0].components observation = model.data[0] - cenx = observation.image.n_cols/2. - ceny = observation.image.n_rows/2. + cenx = observation.image.n_cols / 2.0 + ceny = observation.image.n_rows / 2.0 if limits_x is not None: limits_x.max = float(observation.image.n_cols) if limits_y is not None: limits_y.max = float(observation.image.n_rows) init_component(comp1, cen_x=cenx, cen_y=ceny) - init_component(comp2, cen_x=cenx, cen_y=ceny, - sigma_x=self.sigma_x, sigma_y=self.sigma_y, rho=self.rho) + init_component( + comp2, cen_x=cenx, cen_y=ceny, sigma_x=self.sigma_x, sigma_y=self.sigma_y, rho=self.rho + ) params_free = get_params_uniq(model, fixed=False) for param in params_free: if isinstance(param, g2f.IntegralParameterD): param.value = self.flux - # We should have done this in get_source_observation, but it gets called first + # Should be done in get_source_observation, but it gets called first + # ... and therefore does not have the initialization above model.setup_evaluators(evaluatormode=g2f.Model.EvaluatorMode.image) model.evaluate() - rng = np.random.default_rng(source['id'] + 100000) + rng = np.random.default_rng(source["id"] + 100000) for idx_obs, observation in enumerate(model.data): image, sigma_inv = observation.image, observation.sigma_inv image.data.flat = model.outputs[0].data sigma_inv.data.flat = np.sqrt(image.data + self.background) - image.data.flat = image.data + sigma_inv.data*rng.standard_normal(image.data.shape) - sigma_inv.data.flat = 1/sigma_inv.data + image.data.flat = image.data + sigma_inv.data * rng.standard_normal(image.data.shape) + sigma_inv.data.flat = 1 / sigma_inv.data # This is mandatory because C++ construction does no initialization # (could instead initialize in get_source_observation) # TODO: Do some timings to see which is more efficient diff --git a/multiprofit/fit_catalog.py b/python/lsst/multiprofit/fit_catalog.py similarity index 74% rename from multiprofit/fit_catalog.py rename to python/lsst/multiprofit/fit_catalog.py index 86c3d36..2d3606c 100644 --- a/multiprofit/fit_catalog.py +++ b/python/lsst/multiprofit/fit_catalog.py @@ -27,29 +27,36 @@ from pydantic.dataclasses import dataclass from .modeller import ModelFitConfig +from .utils import ArbitraryAllowedConfig class CatalogExposureABC(ABC): + """Interface for catalog-exposure pairs.""" + @abstractmethod def get_catalog(self) -> Iterable: """Return a row-iterable catalog covering an exposure.""" -class ColumnInfoConfig: - arbitrary_types_allowed = True +# class ColumnInfoConfig: +# """Pydantic config to allow arbitrary typed Fields.""" +# +# arbitrary_types_allowed = True -@dataclass(frozen=True, kw_only=True, config=ColumnInfoConfig) +@dataclass(frozen=True, kw_only=True, config=ArbitraryAllowedConfig) class ColumnInfo: """Metadata for a column in a catalog.""" + dtype: str = pydantic.Field(title="Column data type name (numpy or otherwise)") key: str = pydantic.Field(title="Column key (name)") description: str = pydantic.Field("", title="Column description") - unit: u.UnitBase = pydantic.Field(u.Unit(""), title="Column unit (astropy)") + unit: u.UnitBase | None = pydantic.Field(None, title="Column unit (astropy)") class CatalogFitterConfig(pexConfig.Config): """Configuration for generic MultiProFit fitting tasks.""" + column_id = pexConfig.Field[str](default="id", doc="Catalog index column key") compute_errors = pexConfig.ChoiceField[str]( default="INV_HESSIAN_BESTFIT", @@ -58,7 +65,7 @@ class CatalogFitterConfig(pexConfig.Config): "NONE": "no errors computed", "INV_HESSIAN": "inverse hessian using noisy image as data", "INV_HESSIAN_BESTFIT": "inverse hessian using best-fit model as data", - } + }, ) compute_errors_from_jacobian = pexConfig.Field[bool]( default=True, @@ -72,8 +79,12 @@ class CatalogFitterConfig(pexConfig.Config): fit_centroid = pexConfig.Field[bool](default=True, doc="Fit centroid parameters") fit_linear_init = pexConfig.Field[bool](default=True, doc="Fit linear parameters after initialization") fit_linear_final = pexConfig.Field[bool](default=True, doc="Fit linear parameters after optimization") - flag_errors = pexConfig.DictField(default={}, keytype=str, itemtype=str, - doc="Flag column names to set, keyed by name of exception to catch") + flag_errors = pexConfig.DictField( + default={}, + keytype=str, + itemtype=str, + doc="Flag column names to set, keyed by name of exception to catch", + ) prefix_column = pexConfig.Field[str](default="mpf_", doc="Column name prefix") def schema( @@ -93,19 +104,21 @@ def schema( An ordered list of ColumnInfo instances. """ schema = [ - ColumnInfo(key=self.column_id, dtype='i8'), - ColumnInfo(key='n_iter', dtype='i4'), - ColumnInfo(key='time_eval', dtype='f8', unit=u.s), - ColumnInfo(key='time_fit', dtype='f8', unit=u.s), - ColumnInfo(key='time_full', dtype='f8', unit=u.s), - ColumnInfo(key='chisq_red', dtype='f8'), - ColumnInfo(key='unknown_flag', dtype='bool'), + ColumnInfo(key=self.column_id, dtype="i8"), + ColumnInfo(key="n_iter", dtype="i4"), + ColumnInfo(key="time_eval", dtype="f8", unit=u.s), + ColumnInfo(key="time_fit", dtype="f8", unit=u.s), + ColumnInfo(key="time_full", dtype="f8", unit=u.s), + ColumnInfo(key="chisq_red", dtype="f8"), + ColumnInfo(key="unknown_flag", dtype="bool"), ] - schema.extend([ColumnInfo(key=key, dtype='bool') for key in self.flag_errors.keys()]) + schema.extend([ColumnInfo(key=key, dtype="bool") for key in self.flag_errors.keys()]) # Always have a centroid column, even if not fitting # It may still be useful for reconstruction - schema.extend([ - ColumnInfo(key='cen_x', dtype='f8', unit=u.pix), - ColumnInfo(key='cen_y', dtype='f8', unit=u.pix), - ]) + schema.extend( + [ + ColumnInfo(key="cen_x", dtype="f8", unit=u.pix), + ColumnInfo(key="cen_y", dtype="f8", unit=u.pix), + ] + ) return schema diff --git a/multiprofit/fit_psf.py b/python/lsst/multiprofit/fit_psf.py similarity index 84% rename from multiprofit/fit_psf.py rename to python/lsst/multiprofit/fit_psf.py index 6505f53..68243e0 100644 --- a/multiprofit/fit_psf.py +++ b/python/lsst/multiprofit/fit_psf.py @@ -39,18 +39,40 @@ from .utils import get_params_uniq +class PsfRebuildFitFlagError(RuntimeError): + """RuntimeError for when a PSF can't be rebuilt because the fit failed""" + + class CatalogExposurePsfABC(CatalogExposureABC): + """A CatalogExposure for PSF fitting.""" + @abstractmethod def get_psf_image(self, source: astropy.table.Row | Mapping[str, Any]) -> np.array: - """Return a normalized, centered, odd-sized PSF image array.""" + """Get a PSF image for a specific source. + + Parameters + ---------- + source + The source row/dict. + + Returns + ------- + The image of the PSF. + + Notes + ----- + The PSF image should be normalized, and centered in a 2D array of odd + dimensions on both sides. + """ class CatalogPsfFitterConfig(CatalogFitterConfig): """Configuration for MultiProFit PSF image fitter.""" + gaussians = pexConfig.ConfigDictField( default={ - 'comp1': GaussianConfig(size=ParameterConfig(value_initial=1.5)), - 'comp2': GaussianConfig(size=ParameterConfig(value_initial=3.0)), + "comp1": GaussianConfig(size=ParameterConfig(value_initial=1.5)), + "comp2": GaussianConfig(size=ParameterConfig(value_initial=3.0)), }, doc="Gaussian components", itemtype=GaussianConfig, @@ -76,12 +98,17 @@ def rebuild_psfmodel(self, params: astropy.table.Row | Mapping[str, Any]) -> g2f psfmodel : `g2f.PsfModel` The rebuilt PSF model. """ + # TODO: Improve _flag checking (add a total _flag column) + for flag in (col for col in params.keys() if col.endswith("_flag")): + if params[flag]: + raise PsfRebuildFitFlagError(f"Failed to rebuild PSF; {flag} set") + n_gaussians = len(self.gaussians) idx_gauss_max = n_gaussians - 1 - sigma_xs = [0.]*n_gaussians - sigma_ys = [0.]*n_gaussians - rhos = [0.]*n_gaussians - fracs = [1.]*n_gaussians + sigma_xs = [0.0] * n_gaussians + sigma_ys = [0.0] * n_gaussians + rhos = [0.0] * n_gaussians + fracs = [1.0] * n_gaussians for idx, (name, config) in enumerate(self.gaussians.items()): sigma_xs[idx] = params[f"{self.prefix_column}{name}_sigma_x"] @@ -102,7 +129,7 @@ def schema( if bands is not None: if len(bands) != 1: raise ValueError("CatalogPsfFitter doesn't support multiple bands") - prefix_band = f'{bands[0]}_' + prefix_band = f"{bands[0]}_" schema = super().schema(bands) n_gaussians = len(self.gaussians) idx_gauss_max = n_gaussians - 1 @@ -110,12 +137,12 @@ def schema( for idx_gauss, name in enumerate(self.gaussians.keys()): prefix_comp = f"{name}_{prefix_band}" columns_comp = [ - ColumnInfo(key=f'{prefix_comp}sigma_x', dtype='f8', unit=u.pix), - ColumnInfo(key=f'{prefix_comp}sigma_y', dtype='f8', unit=u.pix), - ColumnInfo(key=f'{prefix_comp}rho', dtype='f8', unit=u.pix), + ColumnInfo(key=f"{prefix_comp}sigma_x", dtype="f8", unit=u.pix), + ColumnInfo(key=f"{prefix_comp}sigma_y", dtype="f8", unit=u.pix), + ColumnInfo(key=f"{prefix_comp}rho", dtype="f8", unit=u.pix), ] if idx_gauss != idx_gauss_max: - columns_comp.append(ColumnInfo(key=f'{prefix_comp}fluxfrac', dtype='f8')) + columns_comp.append(ColumnInfo(key=f"{prefix_comp}fluxfrac", dtype="f8")) schema.extend(columns_comp) return schema @@ -141,6 +168,7 @@ class CatalogPsfFitter: Any exceptions raised and not in errors_expected will be logged in a generic unknown_flag failure column. """ + def __init__( self, modeller: Modeller = None, @@ -172,20 +200,22 @@ def _get_data(img_psf: np.array, gain: float = 1e5) -> g2f.Data: """ # TODO: Improve these arbitrary definitions # Look at S/N of PSF stars? - background = np.std(img_psf[img_psf < 2*np.abs(np.min(img_psf))]) + background = np.std(img_psf[img_psf < 2 * np.abs(np.min(img_psf))]) # Hacky fix; PSFs shouldn't have negative values but often do min_psf = np.min(img_psf) if not (background > -min_psf): background = -1.1 * min_psf img_sig_inv = np.sqrt(gain / (img_psf + background)) - return g2f.Data([ - g2f.Observation( - channel=g2f.Channel.NONE, - image=g2.ImageD(img_psf), - sigma_inv=g2.ImageD(img_sig_inv), - mask_inv=g2.ImageB(np.ones_like(img_psf)), - ) - ]) + return g2f.Data( + [ + g2f.Observation( + channel=g2f.Channel.NONE, + image=g2.ImageD(img_psf), + sigma_inv=g2.ImageD(img_sig_inv), + mask_inv=g2.ImageB(np.ones_like(img_psf)), + ) + ] + ) @staticmethod def _get_logger() -> logging.Logger: @@ -201,7 +231,7 @@ def fit( catexp: CatalogExposurePsfABC, config: CatalogPsfFitterConfig = None, logger: logging.Logger = None, - **kwargs + **kwargs, ) -> astropy.table.Table: """Fit PSF models for a catalog with MultiProFit. @@ -261,7 +291,8 @@ def fit( flux_total = flux_total[0] # TODO: Remove isinstance when channel filtering is fixed fluxfracs = tuple( - x for x in get_params_uniq(model_source, linear=False, channel=g2f.Channel.NONE, fixed=False) + x + for x in get_params_uniq(model_source, linear=False, channel=g2f.Channel.NONE, fixed=False) if isinstance(x, g2f.ProperFractionParameterD) ) if len(fluxfracs) != (n_gaussians - 1): @@ -283,9 +314,10 @@ def fit( idx_var_first = keys.index("cen_x") columns_write = [f"{prefix}{col.key}" for col in columns[idx_var_first:]] dtypes = [(f'{prefix if col.key != config.column_id else ""}{col.key}', col.dtype) for col in columns] - meta = {'config': config.toDict()} - results = Table(data=np.full(n_rows, np.nan, dtype=dtypes), units=[x.unit for x in columns], - meta=meta) + meta = {"config": config.toDict()} + results = Table( + data=np.full(n_rows, np.nan, dtype=dtypes), units=[x.unit for x in columns], meta=meta + ) # Set nan-default flags to False instead for flag in columns[idx_flag_first:idx_var_first]: results[f"{prefix}{flag.key}"] = False @@ -303,17 +335,17 @@ def fit( try: img_psf = catexp.get_psf_image(source) - cenx.value, ceny.value = (x/2. for x in img_psf.shape[::-1]) + cenx.value, ceny.value = (x / 2.0 for x in img_psf.shape[::-1]) + data = CatalogPsfFitter._get_data(img_psf) + model = g2f.Model(data=data, psfmodels=[model_psf], sources=[model_source], priors=priors) + self.initialize_model(model=model, config=config, source=source) + # Caches the jacobian residual if the kernel size is unchanged if img_psf.size != size: fitInputs = None size = np.copy(img_psf.size) else: fitInputs = fitInputs if not fitInputs.validate_for_model(model) else None - - data = CatalogPsfFitter._get_data(img_psf) - model = g2f.Model(data=data, psfmodels=[model_psf], sources=[model_source], priors=priors) - self.initialize_model(model=model, config=config, source=source) if config.fit_linear_init: flux_total.fixed = False @@ -322,20 +354,23 @@ def fit( result = self.modeller.fit_gaussians_linear(gaussians_linear, data[0]) result = list(result.values())[0] # Re-normalize fluxes (hopefully close already) - result = np.clip(result*np.array([ - x[0].at(0).integral.value for x in gaussians_linear.gaussians_free - ]), 1e-2, 0.99) + result = np.clip( + result + * np.array([x[0].at(0).integral.value for x in gaussians_linear.gaussians_free]), + 1e-2, + 0.99, + ) result /= np.sum(result) for idx_param, param in enumerate(fluxfracs): param.value = result[idx_param] - result /= np.sum(result[idx_param + 1:]) + result /= np.sum(result[idx_param + 1 :]) result_full = self.modeller.fit_model(model, fitinputs=fitInputs, **kwargs) fitInputs = result_full.inputs results[f"{prefix}n_iter"][idx] = result_full.n_eval_func results[f"{prefix}time_eval"][idx] = result_full.time_eval results[f"{prefix}time_fit"][idx] = result_full.time_run - results[f"{prefix}chisq_red"][idx] = np.sum(fitInputs.residual**2)/size + results[f"{prefix}chisq_red"][idx] = np.sum(fitInputs.residual**2) / size for param, value, column in zip(params, result_full.params_best, columns_write): param.value_transformed = value @@ -347,8 +382,9 @@ def fit( column = self.errors_expected.get(e.__class__, "") if column: row[f"{prefix}{column}"] = True - logger.info(f"{id_source=} ({idx}/{n_rows}) PSF fit failed with not unexpected" - f" exception={e}") + logger.info( + f"{id_source=} ({idx}/{n_rows}) PSF fit failed with not unexpected" f" exception={e}" + ) else: row[f"{prefix}unknown_flag"] = True logger.info(f"{id_source=} ({idx}/{n_rows}) PSF fit failed with unexpected exception={e}") @@ -378,7 +414,7 @@ def initialize_model( Hard limits for the source's y centroid. """ n_rows, n_cols = model.data[0].image.data.shape - cen_x, cen_y = n_cols/2.0, n_rows/2.0 + cen_x, cen_y = n_cols / 2.0, n_rows / 2.0 centroids = set() if limits_x is None: limits_x = g2f.LimitsD(0, n_cols) diff --git a/multiprofit/fit_source.py b/python/lsst/multiprofit/fit_source.py similarity index 78% rename from multiprofit/fit_source.py rename to python/lsst/multiprofit/fit_source.py index 40ba3d9..b54d5f9 100644 --- a/multiprofit/fit_source.py +++ b/python/lsst/multiprofit/fit_source.py @@ -40,6 +40,8 @@ class CatalogExposureSourcesABC(CatalogExposureABC): + """Interface for a CatalogExposure for source modelling.""" + @property def band(self) -> str: """Return the name of the exposure's passband (e.g. 'r').""" @@ -87,14 +89,17 @@ def get_source_observation(self, source: Mapping[str, Any]) -> g2f.Observation: class CatalogSourceFitterConfig(CatalogFitterConfig): """Configuration for the MultiProFit profile fitter.""" + convert_cen_xy_to_radec = pexConfig.Field[bool](default=True, doc="Convert cen x/y params to RA/dec") fit_cen_x = pexConfig.Field[bool](default=True, doc="Fit x centroid parameter") fit_cen_y = pexConfig.Field[bool](default=True, doc="Fit y centroid parameter") n_pointsources = pexConfig.Field[int](default=0, doc="Number of central point source components") - prior_cen_x_stddev = pexConfig.Field[float](default=0, - doc="Prior std. dev. on x centroid (ignored if not >0)") - prior_cen_y_stddev = pexConfig.Field[float](default=0, - doc="Prior std. dev. on y centroid (ignored if not >0)") + prior_cen_x_stddev = pexConfig.Field[float]( + default=0, doc="Prior std. dev. on x centroid (ignored if not >0)" + ) + prior_cen_y_stddev = pexConfig.Field[float]( + default=0, doc="Prior std. dev. on y centroid (ignored if not >0)" + ) # TODO: Verify that component names don't clash sersics = pexConfig.ConfigDictField( @@ -117,8 +122,8 @@ def make_model_data( else: model_source, priors = model_priors n_catexps = len(catexps) - observations = [None]*n_catexps - psfmodels = [None]*n_catexps + observations = [None] * n_catexps + psfmodels = [None] * n_catexps for idx_catexp in range(n_catexps): catexp = catexps[idx_catexp] @@ -136,8 +141,14 @@ def make_model_data( def make_source( self, channels: list[g2f.Channel], - ) -> tuple[g2f.Source, list[g2f.Prior], g2f.LimitsD, g2f.LimitsD, - g2f.CentroidXParameterD, g2f.CentroidXParameterD]: + ) -> tuple[ + g2f.Source, + list[g2f.Prior], + g2f.LimitsD, + g2f.LimitsD, + g2f.CentroidXParameterD, + g2f.CentroidXParameterD, + ]: limits_x = g2f.LimitsD(min=0, max=np.Inf) limits_y = g2f.LimitsD(min=0, max=np.Inf) cen_x = g2f.CentroidXParameterD(value=0, limits=limits_x, fixed=not self.fit_cen_x) @@ -145,7 +156,7 @@ def make_source( centroid = g2f.CentroidParameters(cen_x, cen_y) priors = [] n_sersics = len(self.sersics) - components = [None]*(self.n_pointsources + n_sersics) + components = [None] * (self.n_pointsources + n_sersics) for idx in range(self.n_pointsources): components[idx] = g2f.GaussianComponent( centroid=centroid, @@ -154,12 +165,17 @@ def make_source( sigma_y=g2f.SigmaYParameterD(0, fixed=True), rho=g2f.RhoParameterD(0, fixed=True), ), - integral=g2f.LinearIntegralModel([ - (channel, g2f.IntegralParameterD( - 1.0, transform=transforms_ref['log10'], label=f'Pt. Src {channel.name}-band' - )) - for channel in channels - ]), + integral=g2f.LinearIntegralModel( + [ + ( + channel, + g2f.IntegralParameterD( + 1.0, transform=transforms_ref["log10"], label=f"Pt. Src {channel.name}-band" + ), + ) + for channel in channels + ] + ), ) idx = self.n_pointsources for sersic in self.sersics.values(): @@ -168,13 +184,9 @@ def make_source( priors.extend(priors_comp) idx += 1 if self.prior_cen_x_stddev > 0 and np.isfinite(self.prior_cen_x_stddev): - priors.append( - g2f.GaussianPrior(centroid.x_param_ptr, 0, self.prior_cen_x_stddev) - ) + priors.append(g2f.GaussianPrior(centroid.x_param_ptr, 0, self.prior_cen_x_stddev)) if self.prior_cen_y_stddev > 0 and np.isfinite(self.prior_cen_y_stddev): - priors.append( - g2f.GaussianPrior(centroid.y_param_ptr, 0, self.prior_cen_y_stddev) - ) + priors.append(g2f.GaussianPrior(centroid.y_param_ptr, 0, self.prior_cen_y_stddev)) return g2f.Source(components), priors, limits_x, limits_y, cen_x, cen_y def schema( @@ -184,40 +196,42 @@ def schema( if bands is None or not (len(bands) > 0): raise ValueError("PSF CatalogSourceFitter must provide at least one band") columns = super().schema(bands) - suffixes = ('', '_err') if self.compute_errors else ('',) + suffixes = ("", "_err") if self.compute_errors else ("",) for suffix in suffixes: - if suffix == '_err': + if suffix == "_err": if self.fit_cen_x: - columns.append(ColumnInfo(key=f'cen_x{suffix}', dtype='f8', unit=u.pix)) + columns.append(ColumnInfo(key=f"cen_x{suffix}", dtype="f8", unit=u.pix)) if self.fit_cen_y: - columns.append(ColumnInfo(key=f'cen_y{suffix}', dtype='f8', unit=u.pix)) + columns.append(ColumnInfo(key=f"cen_y{suffix}", dtype="f8", unit=u.pix)) for idx in range(self.n_pointsources): prefix_comp = f"ps{idx + 1}_" for band in bands: - columns.append(ColumnInfo(key=f'{prefix_comp}{band}_flux{suffix}', dtype='f8')) + columns.append(ColumnInfo(key=f"{prefix_comp}{band}_flux{suffix}", dtype="f8")) for idx, name_comp in enumerate(self.sersics.keys()): prefix_comp = f"{name_comp}_" columns_comp = [ - ColumnInfo(key=f'{prefix_comp}sigma_x{suffix}', dtype='f8', unit=u.pix), - ColumnInfo(key=f'{prefix_comp}sigma_y{suffix}', dtype='f8', unit=u.pix), - ColumnInfo(key=f'{prefix_comp}rho{suffix}', dtype='f8', unit=u.pix), + ColumnInfo(key=f"{prefix_comp}sigma_x{suffix}", dtype="f8", unit=u.pix), + ColumnInfo(key=f"{prefix_comp}sigma_y{suffix}", dtype="f8", unit=u.pix), + ColumnInfo(key=f"{prefix_comp}rho{suffix}", dtype="f8", unit=u.pix), ] for band in bands: columns_comp.append( - ColumnInfo(key=f'{prefix_comp}{band}_flux{suffix}', dtype='f8', - unit=u.Unit(self.unit_flux)) + ColumnInfo( + key=f"{prefix_comp}{band}_flux{suffix}", dtype="f8", + unit=u.Unit(self.unit_flux) if self.unit_flux else None, + ) ) columns.extend(columns_comp) if self.convert_cen_xy_to_radec: - columns.append(ColumnInfo(key=f'cen_ra{suffix}', dtype='f8', unit=u.deg)) - columns.append(ColumnInfo(key=f'cen_dec{suffix}', dtype='f8', unit=u.deg)) + columns.append(ColumnInfo(key=f"cen_ra{suffix}", dtype="f8", unit=u.deg)) + columns.append(ColumnInfo(key=f"cen_dec{suffix}", dtype="f8", unit=u.deg)) return columns @dataclass(kw_only=True) class CatalogSourceFitterABC(ABC): - """ Fit a Gaussian mixture source model to an image with a PSF model. + """Fit a Gaussian mixture source model to an image with a PSF model. Parameters ---------- @@ -251,7 +265,7 @@ def fit( catexps: list[CatalogExposureSourcesABC], config: CatalogSourceFitterConfig = None, logger: logging.Logger = None, - **kwargs + **kwargs, ) -> astropy.table.Table: """Fit PSF-convolved source models with MultiProFit. @@ -315,16 +329,15 @@ def fit( idx_var_last = idx_var_first + len(params) column_cen_x, column_cen_y = (f"{prefix}cen_{v}" for v in ("x", "y")) columns_write = [f"{prefix}{col.key}" for col in columns[idx_var_first:idx_var_last]] - n_radec = 2*config.convert_cen_xy_to_radec - columns_radec = [f"{prefix}{col.key}" for col in columns[idx_var_last:idx_var_last + n_radec]] + n_radec = 2 * config.convert_cen_xy_to_radec + columns_radec = [f"{prefix}{col.key}" for col in columns[idx_var_last : idx_var_last + n_radec]] idx_var_last += n_radec - columns_write_err = [f"{prefix}{col.key}" for col in columns[idx_var_last:len(columns) - n_radec]] + columns_write_err = [f"{prefix}{col.key}" for col in columns[idx_var_last : len(columns) - n_radec]] assert len(columns_write_err) == len(params) - columns_radec_err = [f"{prefix}{col.key}" for col in columns[len(columns) - n_radec:len(columns)]] + columns_radec_err = [f"{prefix}{col.key}" for col in columns[len(columns) - n_radec : len(columns)]] dtypes = [(f'{prefix if col.key != config.column_id else ""}{col.key}', col.dtype) for col in columns] - meta = {'config': config.toDict()} - results = Table(data=np.full(n_rows, 0, dtype=dtypes), units=[x.unit for x in columns], - meta=meta) + meta = {"config": config.toDict()} + results = Table(data=np.full(n_rows, 0, dtype=dtypes), units=[x.unit for x in columns], meta=meta) # Validate that the columns are in the right order # assert because this is a logic error if it fails @@ -347,10 +360,10 @@ def fit( kwargs_err_default = { True: { - 'options': g2f.HessianOptions(findiff_add=1e-3, findiff_frac=1e-3), - 'use_diag_only': config.compute_errors_no_covar, + "options": g2f.HessianOptions(findiff_add=1e-3, findiff_frac=1e-3), + "use_diag_only": config.compute_errors_no_covar, }, - False: {'options': g2f.HessianOptions(findiff_add=1e-6, findiff_frac=1e-6)}, + False: {"options": g2f.HessianOptions(findiff_add=1e-6, findiff_frac=1e-6)}, } for idx in range_idx: @@ -362,7 +375,9 @@ def fit( try: model, data, psfmodels = config.make_model_data( - idx_source=idx, model_priors=(model_source, priors), catexps=catexps, + idx_source=idx, + model_priors=(model_source, priors), + catexps=catexps, ) self.initialize_model(model, source_multi, limits_x=limits_x, limits_y=limits_y) @@ -376,7 +391,8 @@ def fit( else: fitInputs = fitInputs if not fitInputs.validate_for_model(model) else None - # TODO: set limits_flux, transform_flux if not config.fit_linear_init + # TODO: set limits_flux, transform_flux + # (if not config.fit_linear_init) if config.fit_linear_init: self.modeller.fit_model_linear(model=model, ratio_min=1e-3) @@ -388,7 +404,8 @@ def fit( results[f"{prefix}time_eval"][idx] = result_full.time_eval results[f"{prefix}time_fit"][idx] = result_full.time_run - # Set all params to best fit values (in case the optimizer didn't) + # Set all params to best fit values + # In case the optimizer doesn't for param, value, column in zip(params, result_full.params_best, columns_write): param.value_transformed = value results[column][idx] = param.value @@ -397,7 +414,7 @@ def fit( loglike_init, loglike_new = self.modeller.fit_model_linear( model=model, ratio_min=0.01, validate=True ) - loglike_diff = np.sum(loglike_new) - np.sum(loglike_init) + np.sum(loglike_new) - np.sum(loglike_init) # TODO: See if it makes sense to set only flux params for param, column in zip(params, columns_write): results[column][idx] = param.value @@ -418,51 +435,62 @@ def fit( errors_iter = None if config.compute_errors_from_jacobian: try: - errors_iter = np.sqrt(self.modeller.compute_variances( - model_eval, transformed=False, use_diag_only=config.compute_errors_no_covar, - )) + errors_iter = np.sqrt( + self.modeller.compute_variances( + model_eval, + transformed=False, + use_diag_only=config.compute_errors_no_covar, + ) + ) errors.append((errors_iter, np.sum(~(errors_iter > 0)))) - except Exception as e: + except Exception: pass if errors_iter is None: + img_data_old = [] if errors_hessian_bestfit: # Model sans prior - model_eval = g2f.Model(data=model.data, psfmodels=model.psfmodels, - sources=model.sources) + model_eval = g2f.Model( + data=model.data, psfmodels=model.psfmodels, sources=model.sources + ) model_eval.setup_evaluators(evaluatormode=model.EvaluatorMode.image) model_eval.evaluate() - img_data_old = [] for obs, output in zip(model_eval.data, model_eval.outputs): img_data_old.append(obs.image.data.copy()) img = obs.image.data - img.flat = ( - output.data.flat - ) - # To make this a real bootstrap, could do this (but would need to iterate): - # + rng.standard_normal(img.size) * obs.sigma_inv.data.flat + img.flat = output.data.flat + # To make this a real bootstrap, could do this + # (but would need to iterate): + # + rng.standard_normal(img.size)*( + # obs.sigma_inv.data.flat) for return_negative in (False, True): kwargs_err = kwargs_err_default[return_negative] if errors and errors[-1][1] == 0: break try: - errors_iter = np.sqrt(self.modeller.compute_variances( - model_eval, transformed=False, **kwargs_err - )) + errors_iter = np.sqrt( + self.modeller.compute_variances( + model_eval, transformed=False, **kwargs_err + ) + ) errors.append((errors_iter, np.sum(~(errors_iter > 0)))) - except Exception as e: + except Exception: try: - errors_iter = np.sqrt(self.modeller.compute_variances( - model_eval, transformed=False, use_svd=True, **kwargs_err, - )) + errors_iter = np.sqrt( + self.modeller.compute_variances( + model_eval, + transformed=False, + use_svd=True, + **kwargs_err, + ) + ) errors.append((errors_iter, np.sum(~(errors_iter > 0)))) except Exception: pass - - if errors_hessian_bestfit: - for obs, img_datum_old in zip(model.data, img_data_old): - obs.image.data.flat = img_datum_old.flat + if errors_hessian_bestfit: + for obs, img_datum_old in zip(model.data, img_data_old): + obs.image.data.flat = img_datum_old.flat if errors: idx_min = np.argmax([err[1] for err in errors]) @@ -470,10 +498,11 @@ def fit( if plot: errors_plot = np.clip(errors, 0, 1000) errors_plot[~np.isfinite(errors_plot)] = 0 - from multiprofit.plots import plot_loglike, plt + from .plots import plot_loglike, ErrorValues + try: - plot_loglike(model, errors=errors) - except Exception as e: + plot_loglike(model, errors={'err': ErrorValues(values=errors_plot)}) + except Exception: for param in params: param.fixed = False @@ -481,20 +510,26 @@ def fit( results[column][idx] = value if config.convert_cen_xy_to_radec: if not config.fit_cen_x and config.fit_cen_y: - # Note this will make naive tests that check for all-finite results fail + # Note this will make naive tests that check + # for all-finite results fail radec_err = np.nan, np.nan else: # TODO: improve this - # For one, it won't work right at RA = 359.9999... - # Could consider dividing by sqrt(2) but it would multiply out later - radec_err = np.array(self.get_model_radec( - source_multi, cen_x.value + errors[0], cen_y.value + errors[1], - )) + # For one, it won't work right at RA=359.99... + # Could consider dividing by sqrt(2) + # but it would multiply out later + radec_err = np.array( + self.get_model_radec( + source_multi, + cen_x.value + errors[0], + cen_y.value + errors[1], + ) + ) radec_err -= radec for value, column in zip(radec_err, columns_radec_err): results[column][idx] = np.abs(value) - results[f"{prefix}chisq_red"][idx] = np.sum(fitInputs.residual**2)/size + results[f"{prefix}chisq_red"][idx] = np.sum(fitInputs.residual**2) / size results[f"{prefix}time_full"][idx] = time.process_time() - time_init except Exception as e: size = 0 if fitInputs is None else size_new @@ -543,9 +578,11 @@ def get_model( catexps The catalog-exposure pairs to reconstruct the model for. config - The configuration used to generate sources. Default-initialized if None. + The configuration used to generate sources. + Default-initialized if None. results - The corresponding best-fit parameter catalog to initialize parameter values from. + The corresponding best-fit parameter catalog to initialize + parameter values from. Returns ------- @@ -567,7 +604,9 @@ def get_model( source_multi = catalog_multi[idx_row] model, data, psfmodels = config.make_model_data( - idx_source=idx_row, model_priors=(model_source, priors), catexps=catexps, + idx_source=idx_row, + model_priors=(model_source, priors), + catexps=catexps, ) self.initialize_model(model, source_multi, limits_x=limits_x, limits_y=limits_y) @@ -577,7 +616,7 @@ def get_model( row = results[idx_row] idx_col_start = columns.index(f"{config.prefix_column}cen_x") n_params = len(params) - for param, column in zip(params, columns[idx_col_start: idx_col_start + n_params]): + for param, column in zip(params, columns[idx_col_start : idx_col_start + n_params]): param.value = row[column] return model diff --git a/multiprofit/limits.py b/python/lsst/multiprofit/limits.py similarity index 100% rename from multiprofit/limits.py rename to python/lsst/multiprofit/limits.py diff --git a/multiprofit/modeller.py b/python/lsst/multiprofit/modeller.py similarity index 85% rename from multiprofit/modeller.py rename to python/lsst/multiprofit/modeller.py index c9014c7..14e21cf 100644 --- a/multiprofit/modeller.py +++ b/python/lsst/multiprofit/modeller.py @@ -34,32 +34,35 @@ from .utils import ArbitraryAllowedConfig, get_params_uniq try: - import fastnnls + # TODO: try importlib.util.find_spec + import fastnnls # noqa + has_fastnnls = True except ImportError: has_fastnnls = False class InvalidProposalError(ValueError): - pass + """Error for an invalid parameter proposal.""" fitmethods_linear = { - 'scipy.optimize.nnls': {}, - 'scipy.optimize.lsq_linear': {'bounds': (1e-5, np.Inf), 'method': 'bvls'}, - 'numpy.linalg.lstsq': {'rcond': 1e-3}, + "scipy.optimize.nnls": {}, + "scipy.optimize.lsq_linear": {"bounds": (1e-5, np.Inf), "method": "bvls"}, + "numpy.linalg.lstsq": {"rcond": 1e-3}, } if has_fastnnls: - fitmethods_linear['fastnnls.fnnls'] = {} + fitmethods_linear["fastnnls.fnnls"] = {} @dataclass(frozen=True, kw_only=True, config=ArbitraryAllowedConfig) class LinearGaussians: - """Helper for linear least-squares fitting of Gaussian mixtures. - """ + """Helper for linear least-squares fitting of Gaussian mixtures.""" + gaussians_fixed: g2.Gaussians = pydantic.Field(title="Fixed Gaussian components") gaussians_free: tuple[tuple[g2.Gaussians, g2f.ParameterD], ...] = pydantic.Field( - title="Free Gaussian components") + title="Free Gaussian components" + ) @staticmethod def make( @@ -102,9 +105,7 @@ def make( n_g = len(gaussians) if n_g != 1: raise ValueError(f"{component=} has {gaussians=} of len {n_g=}!=1") - param_flux = component.parameters( - paramfilter=g2f.ParamFilter(nonlinear=False, channel=channel) - ) + param_flux = component.parameters(paramfilter=g2f.ParamFilter(nonlinear=False, channel=channel)) if len(param_flux) != 1: raise ValueError(f"Can't make linear source from {component=} with {len(param_flux)=}") param_flux: g2f.ParameterD = param_flux[0] @@ -113,8 +114,9 @@ def make( else: gaussians_free.append((gaussians, param_flux)) - return LinearGaussians(gaussians_fixed=g2.Gaussians(gaussians_fixed), - gaussians_free=tuple(gaussians_free)) + return LinearGaussians( + gaussians_fixed=g2.Gaussians(gaussians_fixed), gaussians_free=tuple(gaussians_free) + ) def make_image_gaussians( @@ -143,11 +145,13 @@ def make_image_gaussians( gaussians_kernel = g2.Gaussians([g2.Gaussian()]) n_gaussians_kernel = len(gaussians_kernel) n_gaussians_source = len(gaussians_source) - gaussians = g2.ConvolvedGaussians([ - g2.ConvolvedGaussian(source=source, kernel=kernel) - for source in (gaussians_source.at(idx) for idx in range(n_gaussians_source)) - for kernel in (gaussians_kernel.at(idx) for idx in range(n_gaussians_kernel)) - ]) + gaussians = g2.ConvolvedGaussians( + [ + g2.ConvolvedGaussian(source=source, kernel=kernel) + for source in (gaussians_source.at(idx) for idx in range(n_gaussians_source)) + for kernel in (gaussians_kernel.at(idx) for idx in range(n_gaussians_kernel)) + ] + ) return g2.make_gaussians_pixel_D(gaussians=gaussians, **kwargs) @@ -164,18 +168,41 @@ def make_psfmodel_null() -> g2f.PsfModel: class FitInputsBase(ABC): + """Interface for inputs to a model fit.""" + @abstractmethod def validate_for_model(self, model: g2f.Model) -> list[str]: - """Validate that existing FitInputs are valid for a Model""" + """Check that this FitInputs is valid for a Model. + + Parameters + ---------- + model + The model to validate with. + + Returns + ------- + errors + A list of validation errors, if any. + """ class FitInputsDummy(FitInputsBase): + """A dummy FitInputs that always fails to validate. + + This class can be used to initialize a FitInputsBase that may be + reassigned to a non-dummy derived instance in a loop. + """ + def validate_for_model(self, model: g2f.Model) -> list[str]: - return ["This is a dummy FitInputs and will never validate",] + return [ + "This is a dummy FitInputs and will never validate", + ] @dataclass(kw_only=True, config=ArbitraryAllowedConfig) class FitInputs(FitInputsBase): + """Model fit inputs for gauss2dfit.""" + jacobian: np.ndarray = pydantic.Field(None, title="The full Jacobian array") jacobians: list[list[g2.ImageD]] = pydantic.Field( title="Jacobian arrays (views) for each observation", @@ -250,21 +277,21 @@ def from_model( size_data = n_data + n_prior_residuals shape_jacobian = (size_data, n_params_jac) jacobian = np.zeros(shape_jacobian) - jacobians = [None]*n_obs - outputs_prior = [None]*n_params_jac + jacobians = [None] * n_obs + outputs_prior = [None] * n_params_jac for idx in range(n_params_jac): outputs_prior[idx] = g2.ImageD(jacobian[n_data:, idx].view().reshape((1, n_prior_residuals))) residual = np.zeros(size_data) - residuals = [None]*n_obs + residuals = [None] * n_obs residuals_prior = g2.ImageD(residual[n_data:].reshape(1, n_prior_residuals)) offset = 0 for idx_obs in range(n_obs): shape = shapes[idx_obs, :] - size_obs = shape[0]*shape[1] + size_obs = shape[0] * shape[1] end = offset + size_obs - jacobians_obs = [None]*n_params_jac + jacobians_obs = [None] * n_params_jac for idx_jac in range(n_params_jac): jacobians_obs[idx_jac] = g2.ImageD(jacobian[offset:end, idx_jac].view().reshape(shape)) jacobians[idx_obs] = jacobians_obs @@ -273,9 +300,12 @@ def from_model( if offset != n_data_obs[idx_obs]: raise RuntimeError(f"Assigned {offset=} data points != {n_data_obs[idx_obs]=}") return cls( - jacobian=jacobian, jacobians=jacobians, - residual=residual, residuals=residuals, - outputs_prior=outputs_prior, residuals_prior=residuals_prior, + jacobian=jacobian, + jacobians=jacobians, + residual=residual, + residuals=residuals, + outputs_prior=outputs_prior, + residuals_prior=residuals_prior, ) def validate_for_model(self, model: g2f.Model) -> list[str]: @@ -326,9 +356,11 @@ def validate_for_model(self, model: g2f.Model) -> list[str]: class ModelFitConfig(pexConfig.Config): + """Configuration for model fitting.""" + fit_linear_iter = pexConfig.Field[int]( doc="The number of iterations to wait before performing a linear fit during optimization." - " Default 0 disables the feature.", + " Default 0 disables the feature.", default=0, ) @@ -339,9 +371,9 @@ def validate(self): @dataclass(kw_only=True, config=ArbitraryAllowedConfig) class FitResult: - """Results from a Modeller fit, including metadata. - """ - # TODO: Why does setting this default to ModelFitConfig() cause a circular import?? + """Results from a Modeller fit, including metadata.""" + + # TODO: Why does setting default=ModelFitConfig() cause a circular import? config: ModelFitConfig = pydantic.Field(None, title="The configuration for fitting") inputs: FitInputs | None = pydantic.Field(None, title="The fit input arrays") result: Any | None = pydantic.Field( @@ -352,12 +384,13 @@ class FitResult: None, title="The best-fit parameter array (un-transformed)", ) - n_eval_resid: int = pydantic.Field( - 0, title="Total number of self-reported residual function evaluations") + n_eval_resid: int = pydantic.Field(0, title="Total number of self-reported residual function evaluations") n_eval_func: int = pydantic.Field( - 0, title="Total number of optimizer-reported fitness function evaluations") + 0, title="Total number of optimizer-reported fitness function evaluations" + ) n_eval_jac: int = pydantic.Field( - 0, title="Total number of optimizer-reported Jacobian function evaluations") + 0, title="Total number of optimizer-reported Jacobian function evaluations" + ) time_eval: float = pydantic.Field(0, title="Total runtime spent in model/Jacobian evaluation") time_run: float = pydantic.Field(0, title="Total runtime spent in fitting, excluding initial setup") @@ -370,6 +403,7 @@ class Modeller: logger : `logging.Logger` The logger. Defaults to calling `_getlogger`. """ + def __init__(self, logger=None): if logger is None: logger = self._get_logger() @@ -384,15 +418,10 @@ def _get_logger(): return logger @staticmethod - def compute_variances( - model: g2f.Model, - use_diag_only: bool = False, - use_svd: bool = False, - **kwargs - ): + def compute_variances(model: g2f.Model, use_diag_only: bool = False, use_svd: bool = False, **kwargs): hessian = model.compute_hessian(**kwargs).data if use_diag_only: - return -1/np.diag(hessian) + return -1 / np.diag(hessian) if use_svd: u, s, v = np.linalg.svd(-hessian) inverse = np.dot(v.transpose(), np.dot(np.diag(s**-1), u.transpose())) @@ -432,7 +461,7 @@ def fit_gaussians_linear( if psfmodel is None: psfmodel = make_psfmodel_null() if fitmethods is None: - fitmethods = {'scipy.optimize.nnls': fitmethods_linear['scipy.optimize.nnls']} + fitmethods = {"scipy.optimize.nnls": fitmethods_linear["scipy.optimize.nnls"]} else: for fitmethod in fitmethods: if fitmethod not in fitmethods_linear: @@ -456,7 +485,8 @@ def fit_gaussians_linear( image_fixed = make_image_gaussians( gaussians_source=gaussians_linear.gaussians_fixed, gaussians_kernel=gaussians_psf, - n_rows=shape[0], n_cols=shape[1], + n_rows=shape[0], + n_cols=shape[1], ).data if mask_inv is not None: image_fixed = image_fixed[mask_inv] @@ -465,22 +495,22 @@ def fit_gaussians_linear( x = np.zeros((size, n_params)) - params = [None]*n_params + params = [None] * n_params for idx_param, (gaussians_free, param) in enumerate(gaussians_linear.gaussians_free): image_free = make_image_gaussians( gaussians_source=gaussians_free, gaussians_kernel=gaussians_psf, - n_rows=shape[0], n_cols=shape[1], + n_rows=shape[0], + n_cols=shape[1], ).data - x[:, idx_param] = ( - (image_free if mask_inv is None else image_free[mask_inv])*sigma_inv - ).flat + x[:, idx_param] = ((image_free if mask_inv is None else image_free[mask_inv]) * sigma_inv).flat params[idx_param] = param y = observation.image.data if plot: import matplotlib.pyplot as plt - plt.imshow(y, origin='lower') + + plt.imshow(y, origin="lower") plt.show() if mask_inv is not None: y = y[mask_inv] @@ -491,15 +521,16 @@ def fit_gaussians_linear( results = {} for fitmethod, kwargs in fitmethods.items(): - if fitmethod == 'scipy.optimize.nnls': + if fitmethod == "scipy.optimize.nnls": values = spopt.nnls(x, y)[0] - elif fitmethod == 'scipy.optimize.lsq_linear': + elif fitmethod == "scipy.optimize.lsq_linear": kwargs = kwargs if kwargs is not None else fitmethods_linear[fitmethod] values = spopt.lsq_linear(x, y, **kwargs).x - elif fitmethod == 'numpy.linalg.lstsq': + elif fitmethod == "numpy.linalg.lstsq": values = np.linalg.lstsq(x, y, **kwargs)[0] - elif fitmethod == 'fastnnls.fnnls': + elif fitmethod == "fastnnls.fnnls": from fastnnls import fnnls + y = x.T.dot(y) x = x.T.dot(x) values = fnnls(x, y) @@ -514,7 +545,7 @@ def fit_model( fitinputs: FitInputs | None = None, printout: bool = False, config: ModelFitConfig = None, - **kwargs + **kwargs, ) -> FitResult: """Fit a model with a nonlinear optimizer. @@ -566,7 +597,7 @@ def residual_func(params_new, model, params, result): if (fit_linear_iter > 0) and ((result.n_eval_resid + 1) % fit_linear_iter == 0): self.fit_model_linear(model, ratio_min=1e-6) time_init = time.process_time() - loglikes = model.evaluate() + model.evaluate() result.time_eval += time.process_time() - time_init result.n_eval_resid += 1 return result.inputs.residual @@ -586,32 +617,38 @@ def jacobian_func(params_new, model, params, result): params_free = tuple(get_params_uniq(model, fixed=False)) n_params_free = len(params_free) - bounds = ([None]*n_params_free, [None]*n_params_free) - params_init = [None]*n_params_free + bounds = ([None] * n_params_free, [None] * n_params_free) + params_init = [None] * n_params_free for idx, param in enumerate(params_free): limits = param.limits - # If the transform has limits and is more restrictive, use those limits - if hasattr(param.transform, 'limits'): + # If the transform has more restrictive limits, use those + if hasattr(param.transform, "limits"): limits_transform = param.transform.limits n_within = limits.check(limits_transform.min) + limits.check(limits_transform.min) if n_within == 2: limits = limits_transform elif n_within != 0: - raise ValueError(f"{param=} {param.limits=} and {param.transform.limits=}" - f" intersect; one must be a subset of the other") + raise ValueError( + f"{param=} {param.limits=} and {param.transform.limits=}" + f" intersect; one must be a subset of the other" + ) bounds[0][idx] = param.transform.forward(limits.min) bounds[1][idx] = param.transform.forward(limits.max) if not limits.check(param.value): - raise RuntimeError(f'{param=}.value_transformed={param.value} not within {limits=}') + raise RuntimeError(f"{param=}.value_transformed={param.value} not within {limits=}") params_init[idx] = param.value_transformed results = FitResult(inputs=fitinputs, config=config) time_init = time.process_time() result_opt = spopt.least_squares( - residual_func, params_init, jac=jacobian_func, bounds=bounds, - args=(model, params_free, results), x_scale='jac', - **kwargs + residual_func, + params_init, + jac=jacobian_func, + bounds=bounds, + args=(model, params_free, results), + x_scale="jac", + **kwargs, ) results.time_run = time.process_time() - time_init results.result = result_opt @@ -658,7 +695,7 @@ def fit_model_linear( values_init[parameter] = float(parameter.value) if not (ratio >= ratio_min): ratio = ratio_min - value_new = max(ratio*parameter.value, parameter.limits.min) + value_new = max(ratio * parameter.value, parameter.limits.min) values_new[parameter] = value_new for parameter, value in values_new.items(): @@ -713,9 +750,11 @@ def make_components_linear( g2f.CentroidXParameterD(gaussian.centroid.x, fixed=True), g2f.CentroidYParameterD(gaussian.centroid.y, fixed=True), ), - g2f.LinearIntegralModel([ - (g2f.Channel.NONE, g2f.IntegralParameterD(gaussian.integral.value)), - ]), + g2f.LinearIntegralModel( + [ + (g2f.Channel.NONE, g2f.IntegralParameterD(gaussian.integral.value)), + ] + ), ) components_new[idx] = component_new return components_new diff --git a/multiprofit/plots.py b/python/lsst/multiprofit/plots.py similarity index 82% rename from multiprofit/plots.py rename to python/lsst/multiprofit/plots.py index f7fddaa..037d70f 100644 --- a/multiprofit/plots.py +++ b/python/lsst/multiprofit/plots.py @@ -14,6 +14,8 @@ @dataclass(kw_only=True) class ErrorValues: + """Configuration for plotting uncertainties.""" + kwargs_plot: dict = field(default_factory=dict) values: np.ndarray @@ -24,7 +26,7 @@ def plot_catalog_bootstrap( paramvals_ref=None, plot_total_fluxes: bool = False, plot_colors: bool = False, - **kwargs + **kwargs, ): """Plot a bootstrap catalog for a single source model. @@ -54,7 +56,7 @@ def plot_catalog_bootstrap( if n_bins is None: n_bins = np.max([int(np.ceil(np.sqrt(n_sources))), 10]) - config = catalog_bootstrap.meta['config'] + config = catalog_bootstrap.meta["config"] prefix = config["prefix_column"] suffix_err = "_err" @@ -65,7 +67,7 @@ def plot_catalog_bootstrap( if paramvals_ref and (len(paramvals_ref) != n_params_init): raise ValueError(f"{len(paramvals_ref)=} != {n_params_init=}") - results_good = catalog_bootstrap[catalog_bootstrap['mpf_n_iter'] > 0] + results_good = catalog_bootstrap[catalog_bootstrap["mpf_n_iter"] > 0] if plot_total_fluxes or plot_colors: if paramvals_ref: @@ -77,14 +79,14 @@ def plot_catalog_bootstrap( results_dict[colname_meas] = results_good[colname_meas] results_dict[colname_err] = results_good[colname_err] - colnames_flux = [colname for colname in colnames_meas if colname.endswith('_flux')] + colnames_flux = [colname for colname in colnames_meas if colname.endswith("_flux")] colnames_flux_band = defaultdict(list) colnames_flux_comp = defaultdict(list) for colname in colnames_flux: colname_short = colname.partition(prefix)[-1] - comp, band = colname_short[:-5].split('_') + comp, band = colname_short[:-5].split("_") colnames_flux_band[band].append(colname) colnames_flux_comp[comp].append(colname) @@ -94,7 +96,7 @@ def plot_catalog_bootstrap( is_err = suffix == "_err" colname_flux = f"{prefix}{band}_flux{suffix}" total = np.sum( - [results_good[f'{colname}{suffix}']**(1 + is_err) for colname in colnames_band], axis=0 + [results_good[f"{colname}{suffix}"] ** (1 + is_err) for colname in colnames_band], axis=0 ) if is_err: total = np.sqrt(total) @@ -106,9 +108,9 @@ def plot_catalog_bootstrap( if band_prev: flux_prev, flux = (results_dict[f"{prefix}{b}_flux"] for b in (band_prev, band)) - mag_prev, mag = (-2.5*np.log10(flux_b) for flux_b in (flux_prev, flux)) + mag_prev, mag = (-2.5 * np.log10(flux_b) for flux_b in (flux_prev, flux)) mag_err_prev, mag_err = ( - results_dict[f"{prefix}{b}_flux{suffix_err}"]/(-0.4*flux_b*ln10) + results_dict[f"{prefix}{b}_flux{suffix_err}"] / (-0.4 * flux_b * ln10) for b, flux_b in ((band_prev, flux_prev), (band, flux)) ) colname_color = f"{prefix}{band_prev}-{band}_flux" @@ -116,10 +118,10 @@ def plot_catalog_bootstrap( colnames_err.append(f"{colname_color}{suffix_err}") results_dict[colname_color] = mag_prev - mag - results_dict[f"{colname_color}{suffix_err}"] = 2.5/ln10*np.hypot(mag_err, mag_err_prev) + results_dict[f"{colname_color}{suffix_err}"] = 2.5 / ln10 * np.hypot(mag_err, mag_err_prev) if paramvals_ref: mag_prev_ref, mag_ref = ( - -2.5*np.log10(paramvals_ref[f"{prefix}{b}_flux"]) for b in (band_prev, band) + -2.5 * np.log10(paramvals_ref[f"{prefix}{b}_flux"]) for b in (band_prev, band) ) paramvals_ref[colname_color] = mag_prev_ref - mag_ref @@ -147,21 +149,21 @@ def plot_catalog_bootstrap( median_err = np.median(errors) axis = ax[idx_row][idx_col] - axis.hist(values, bins=n_bins, color='b', label='fit values', **kwargs) + axis.hist(values, bins=n_bins, color="b", label="fit values", **kwargs) - label = 'median +/- stddev' + label = "median +/- stddev" for offset in (-std, 0, std): - axis.axvline(median + offset, label=label, color='k') + axis.axvline(median + offset, label=label, color="k") label = None if paramvals_ref is not None: value_ref = paramvals_ref[idx_colname] - label_value = f' {value_ref=:.3e} bias={median - value_ref:.3e}' - axis.axvline(value_ref, label='reference', color='k', linestyle='--') + label_value = f" {value_ref=:.3e} bias={median - value_ref:.3e}" + axis.axvline(value_ref, label="reference", color="k", linestyle="--") else: - label_value = f' {median=:.3e}' - axis.hist(median + errors, bins=n_bins, color='r', label='median + error', **kwargs) - axis.set_title(f'{colname_short} {std=:.3e} vs {median_err=:.3e}') - axis.set_xlabel(f'{colname_short} {label_value}') + label_value = f" {median=:.3e}" + axis.hist(median + errors, bins=n_bins, color="r", label="median + error", **kwargs) + axis.set_title(f"{colname_short} {std=:.3e} vs {median_err=:.3e}") + axis.set_xlabel(f"{colname_short} {label_value}") axis.legend() idx_col += 1 @@ -174,8 +176,11 @@ def plot_catalog_bootstrap( def plot_loglike( - model: g2f.Model, params: list[g2f.ParameterD] = None, n_values: int = 15, - errors: dict[str, ErrorValues] = None, values_reference: np.ndarray = None, + model: g2f.Model, + params: list[g2f.ParameterD] = None, + n_values: int = 15, + errors: dict[str, ErrorValues] = None, + values_reference: np.ndarray = None, ): """Plot the loglikehood and derivatives vs free parameter values around best-fit values. @@ -214,12 +219,12 @@ def plot_loglike( raise ValueError(f"{len(values_reference)=} != {n_params=}") n_rows = n_params - fig, ax = plt.subplots(nrows=n_rows, ncols=2, figsize=(10, 3*n_rows)) + fig, ax = plt.subplots(nrows=n_rows, ncols=2, figsize=(10, 3 * n_rows)) axes = [ax] if (n_rows == 1) else ax n_loglikes = len(loglike_init) labels = [channel.name for channel in model.data.channels] - labels.extend(['prior', 'total']) + labels.extend(["prior", "total"]) for param in params: param.fixed = True @@ -228,10 +233,10 @@ def plot_loglike( value_init = param.value param.fixed = False values = [value_init] - loglikes = [loglike_init*0] + loglikes = [loglike_init * 0] dlls = [loglike_grads[row]] - diff_init = 1e-4*np.sign(loglike_grads[row]) + diff_init = 1e-4 * np.sign(loglike_grads[row]) diff = diff_init # TODO: This entire scheme should be improved/replaced @@ -239,7 +244,7 @@ def plot_loglike( # Option: Try to fit a curve once there are a couple of points # on each side of the peak idx_prev = -1 - for idx in range(2*n_values): + for idx in range(2 * n_values): try: param.value_transformed += diff loglikes_new = np.array(model.evaluate()) - loglike_init @@ -258,7 +263,7 @@ def plot_loglike( idx_prev = 0 else: idx_prev = -1 - except RuntimeError as e: + except RuntimeError: break param.value = value_init param.fixed = True @@ -272,28 +277,28 @@ def plot_loglike( for idx in range(n_loglikes): subplot.plot(values, [loglike[idx] for loglike in loglikes], label=labels[idx]) subplot.plot(values, np.sum(loglikes, axis=1), label=labels[-1]) - vline_kwargs = dict(ymin=np.min(loglikes) - 1, ymax=np.max(loglikes) + 1, color='k') + vline_kwargs = dict(ymin=np.min(loglikes) - 1, ymax=np.max(loglikes) + 1, color="k") subplot.vlines(value_init, **vline_kwargs) - suffix = f' {param.label}' if param.label else '' + suffix = f" {param.label}" if param.label else "" subplot.legend() subplot.set_title(f"{param.name}{suffix}") - subplot.set_ylabel('loglike') - subplot.set_ylim(vline_kwargs['ymin'], vline_kwargs['ymax']) + subplot.set_ylabel("loglike") + subplot.set_ylim(vline_kwargs["ymin"], vline_kwargs["ymax"]) subplot = axes[row][1] subplot.plot(values, dlls) - subplot.axhline(0, color='k') + subplot.axhline(0, color="k") subplot.set_ylabel("dloglike/dx") vline_kwargs = dict(ymin=np.min(dlls), ymax=np.max(dlls)) - subplot.vlines(value_init, **vline_kwargs, color='k', label='fit') + subplot.vlines(value_init, **vline_kwargs, color="k", label="fit") if values_reference is not None: - subplot.vlines(values_reference[row], **vline_kwargs, color='b', label='ref') + subplot.vlines(values_reference[row], **vline_kwargs, color="b", label="ref") cycler_linestyle = cycle(linestyles_default) for name_error, valerr in errors.items(): - linestyle = valerr.kwargs_plot.pop('linestyle', next(cycler_linestyle)) + linestyle = valerr.kwargs_plot.pop("linestyle", next(cycler_linestyle)) for idx_ax in range(2): axes[row][idx_ax].vlines( [value_init - valerr.values[row], value_init + valerr.values[row]], diff --git a/multiprofit/priors.py b/python/lsst/multiprofit/priors.py similarity index 64% rename from multiprofit/priors.py rename to python/lsst/multiprofit/priors.py index 5c20176..89a0047 100644 --- a/multiprofit/priors.py +++ b/python/lsst/multiprofit/priors.py @@ -27,6 +27,8 @@ class ShapePriorConfig(pexConfig.Config): + """Configuration for a shape prior.""" + prior_axrat_mean = pexConfig.Field[float]( default=0.7, doc="Prior mean on axis ratio (prior ignored if not >0)", @@ -49,24 +51,40 @@ def get_shape_prior(self, ellipse: g2f.ParametricEllipse) -> g2f.ShapePrior | No use_prior_size = self.prior_size_stddev > 0 and np.isfinite(self.prior_size_stddev) if use_prior_axrat or use_prior_size: - prior_size = g2f.ParametricGaussian1D( - g2f.MeanParameterD(self.prior_size_mean, transform=transforms_ref["log10"]), - g2f.StdDevParameterD(self.prior_size_stddev) - ) if use_prior_size else None - prior_axrat = g2f.ParametricGaussian1D( - g2f.MeanParameterD(self.prior_axrat_mean, transform=transforms_ref["logit_axrat_prior"]), - g2f.StdDevParameterD(self.prior_axrat_stddev) - ) if use_prior_axrat else None + prior_size = ( + g2f.ParametricGaussian1D( + g2f.MeanParameterD(self.prior_size_mean, transform=transforms_ref["log10"]), + g2f.StdDevParameterD(self.prior_size_stddev), + ) + if use_prior_size + else None + ) + prior_axrat = ( + g2f.ParametricGaussian1D( + g2f.MeanParameterD(self.prior_axrat_mean, transform=transforms_ref["logit_axrat_prior"]), + g2f.StdDevParameterD(self.prior_axrat_stddev), + ) + if use_prior_axrat + else None + ) return g2f.ShapePrior(ellipse, prior_size, prior_axrat) return None def get_hst_size_prior(mag_psf_i): - """ Return the mean and stddev for a reasonable HST-based size (major axis half-light radius) prior. + """Return the mean and stddev for an HST-based size prior. + + The size is major axis half-light radius. + + Parameters + ---------- + mag_psf_i + i-band PSF magnitudes of the source(s). Notes ----- - Return values are log10 scaled in units of arcseconds. The input should be a PSF mag because other - magnitudes - even Gaussian - can be unreliable for low S/N (non-)detections. + Return values are log10 scaled in units of arcseconds. + The input should be a PSF mag because other magnitudes - even Gaussian - + can be unreliable for low S/N (non-)detections. """ - return 0.75*(19 - np.clip(mag_psf_i, 10, 30))/6.5, 0.2 + return 0.75 * (19 - np.clip(mag_psf_i, 10, 30)) / 6.5, 0.2 diff --git a/multiprofit/psfmodel_utils.py b/python/lsst/multiprofit/psfmodel_utils.py similarity index 84% rename from multiprofit/psfmodel_utils.py rename to python/lsst/multiprofit/psfmodel_utils.py index 51cdf08..8bf2953 100644 --- a/multiprofit/psfmodel_utils.py +++ b/python/lsst/multiprofit/psfmodel_utils.py @@ -23,7 +23,7 @@ import numpy as np -from multiprofit.limits import limits_ref +from .limits import limits_ref def make_psf_source( @@ -62,7 +62,7 @@ def make_psf_source( Parameter lists must all be the same length. """ if limits_rho is None: - limits_rho = limits_ref['rho'] + limits_rho = limits_ref["rho"] if sigma_xs is None: sigma_xs = [1.5, 3.0] if sigma_ys is not None else sigma_ys if sigma_ys is None: @@ -71,15 +71,15 @@ def make_psf_source( if n_gaussians == 0: raise ValueError(f"{n_gaussians=}!>0") if rhos is None: - rhos = [0.]*n_gaussians + rhos = [0.0] * n_gaussians if fracs is None: - fracs = np.arange(1, n_gaussians + 1)/n_gaussians + fracs = np.arange(1, n_gaussians + 1) / n_gaussians if transforms is None: transforms = {} transforms_default = { - 'frac': transforms.get('frac', g2f.LogitTransformD()), - 'rho': transforms.get('rho', g2f.LogitLimitedTransformD(limits=limits_rho)), - 'sigma': transforms.get('sigma', g2f.Log10TransformD()), + "frac": transforms.get("frac", g2f.LogitTransformD()), + "rho": transforms.get("rho", g2f.LogitLimitedTransformD(limits=limits_rho)), + "sigma": transforms.get("sigma", g2f.Log10TransformD()), } for key, value in transforms_default.items(): if key not in transforms: @@ -98,7 +98,7 @@ def make_psf_source( errors.append(f"fluxes[{idx}]={frac} !>0") if errors: raise ValueError("; ".join(errors)) - fracs[-1] = 1. + fracs[-1] = 1.0 components = [None] * n_gaussians cenx = g2f.CentroidXParameterD(0, limits=g2f.LimitsD(min=0, max=100)) @@ -112,20 +112,21 @@ def make_psf_source( is_last = c == n_last last = g2f.FractionalIntegralModel( [ - (g2f.Channel.NONE, g2f.ProperFractionParameterD( - fracs[c], fixed=is_last, transform=transforms['frac'] - )) + ( + g2f.Channel.NONE, + g2f.ProperFractionParameterD(fracs[c], fixed=is_last, transform=transforms["frac"]), + ) ], - g2f.LinearIntegralModel([ - (g2f.Channel.NONE, g2f.IntegralParameterD(1.0, fixed=True)) - ]) if (c == 0) else last, + g2f.LinearIntegralModel([(g2f.Channel.NONE, g2f.IntegralParameterD(1.0, fixed=True))]) + if (c == 0) + else last, is_last, ) components[c] = g2f.GaussianComponent( g2f.GaussianParametricEllipse( - g2f.SigmaXParameterD(sigma_xs[c], transform=transforms['sigma']), - g2f.SigmaYParameterD(sigma_ys[c], transform=transforms['sigma']), - g2f.RhoParameterD(rhos[c], transform=transforms['rho']), + g2f.SigmaXParameterD(sigma_xs[c], transform=transforms["sigma"]), + g2f.SigmaYParameterD(sigma_ys[c], transform=transforms["sigma"]), + g2f.RhoParameterD(rhos[c], transform=transforms["rho"]), ), centroid, last, diff --git a/multiprofit/transforms.py b/python/lsst/multiprofit/transforms.py similarity index 57% rename from multiprofit/transforms.py rename to python/lsst/multiprofit/transforms.py index 8ce97d5..d845b1a 100644 --- a/multiprofit/transforms.py +++ b/python/lsst/multiprofit/transforms.py @@ -25,11 +25,31 @@ from .limits import limits_ref -def get_logit_limited(lower, upper, factor=1.0, name=None): +def get_logit_limited(lower: float, upper: float, factor: float = 1.0, name: str | None = None): + """Get a logit transform stretched to span a different range than [0,1]. + + Parameters + ---------- + lower + The lower limit of the range to span. + upper + The upper limit of the range to span. + factor + A multiplicative factor to apply to the transformed result. + name + A descriptive name for the transform. + + Returns + ------- + A modified logit transform as specified. + """ return g2f.LogitLimitedTransformD( limits=g2f.LimitsD( - min=lower, max=upper, name=name if name is not None else - f"LogitLimitedTransformD(min={lower}, max={upper}, factor={factor})", + min=lower, + max=upper, + name=name + if name is not None + else f"LogitLimitedTransformD(min={lower}, max={upper}, factor={factor})", ), factor=factor, ) @@ -38,25 +58,44 @@ def get_logit_limited(lower, upper, factor=1.0, name=None): def verify_transform_derivative( transform: g2f.TransformD, value_transformed: float, - derivative: float = None, + derivative: float | None = None, abs_max: float = 1e6, dx_ratios=None, - **kwargs + **kwargs, ): - """ Verify a transform derivative at a given value. + """Verify that the derivative of a transform class is correct. + + Parameters + ---------- + transform + The transform to verify. + value_transformed + The un-transformed value at which to verify the transform. + derivative + The nominal derivative at value_transformed. + Must equal transform.derivative(value_transformed). + abs_max + x value where verification is skipped if np.abs(derivative) > x. + dx_ratios + Iterable of signed ratios to set dx for finite differencing. + dx = value*ratio (untransformed). Only used if dx is None. + kwargs + Keyword arguments to pass to np.isclose when comparing derivatives to + finite differences. + + Raises + ------ + RuntimeError + Raised if the transform derivative doesn't match finite differences + within the specified tolerances. + + Notes + ----- + derivative should only be specified if it has previously been computed for + the exact value_transformed, to avoid re-computing it unnecessarily. - :param transform: gauss2d.fit.TransformD; a transform to verify. - :param value_transformed: float; the un-transformed value at which to verify the transform. - :param abs_max: float; x value where verification skipped if np.abs(derivative) > x - :param dx_ratios: float[]; iterable of ratios to set dx for finite differencing (signed), - where dx = value*ratio (not transformed). Only used if dx is None. - Default: [1e-4, -1e-4, 1e-6, 1e-6, 1e-8, -1e-8, 1e-10, -1e-10, 1e-12, -1e-12, 1e-14, -1e-14] - Verification will test all values in order until at least one passes. - :param kwargs: dict; args to pass to np.isclose when comparing derivative to finite difference e.g. - atol, rtol - :param derivative: float; the nominal derivative at value_transformed. - Must equal transform.derivative(value_transformed) for correct results. - :return: bool; whether the derivative is correct or not. + Default dx_ratios are [1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14]. + Verification will test all ratios until at least one passes. """ # Skip testing finite differencing if the derivative is very large # This might happen e.g. near the limits of the transformation @@ -80,8 +119,9 @@ def verify_transform_derivative( break if not is_close: raise RuntimeError( - f'{transform} derivative={derivative:.8e} != last ' - f'finite diff.={fin_diff:8e} with dx={dx} dx_abs_max={abs_max}') + f"{transform} derivative={derivative:.8e} != last " + f"finite diff.={fin_diff:8e} with dx={dx} dx_abs_max={abs_max}" + ) transforms_ref = { diff --git a/multiprofit/utils.py b/python/lsst/multiprofit/utils.py similarity index 85% rename from multiprofit/utils.py rename to python/lsst/multiprofit/utils.py index bdcbad0..1ea94a3 100644 --- a/multiprofit/utils.py +++ b/python/lsst/multiprofit/utils.py @@ -20,12 +20,18 @@ # along with this program. If not, see . import gauss2d.fit as g2f +import numpy import numpy as np class ArbitraryAllowedConfig: + """Pydantic config to allow arbitrary typed Fields. + + Also disallows unused init kwargs. + """ + arbitrary_types_allowed = True - extra = 'forbid' + extra = "forbid" def get_params_uniq(parametric: g2f.Parametric, **kwargs): @@ -33,7 +39,8 @@ def get_params_uniq(parametric: g2f.Parametric, **kwargs): return {p: None for p in parametric.parameters(paramfilter=g2f.ParamFilter(**kwargs))}.keys() -def normalize(ndarray, return_sum=False): +def normalize(ndarray: numpy.ndarray, return_sum: bool = False): + """Normalize a numpy array.""" total = np.sum(ndarray) ndarray /= total if return_sum: diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 037d3d7..0000000 --- a/setup.cfg +++ /dev/null @@ -1,15 +0,0 @@ -[flake8] -max-line-length = 110 -max-doc-length = 79 -ignore = E133, E226, E228, N802, N803, N806, N812, N813, N815, N816, W503 -exclude = - bin, - doc/conf.py, - **/*/__init__.py, - **/*/version.py, - tests/.tests - -# [tool:pytest] -# addopts = --flake8 -# flake8-ignore = E133 E226 E228 N802 N803 N806 N812 N813 N815 N816 W503 - diff --git a/setup.py b/setup.py deleted file mode 100644 index 269647e..0000000 --- a/setup.py +++ /dev/null @@ -1,42 +0,0 @@ -# This file is part of multiprofit. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import os -import re -import sys -import platform -import subprocess - -from distutils.version import LooseVersion -from setuptools import setup, Extension -from setuptools.command.build_ext import build_ext - -setup( - name='multiprofit', - version='0.1', - author='Dan Taranu', - author_email='dan.s.taranu@gmail.com', - description='Multi-(Object/Band) Source Profile Fitting for Astronomy', - long_description='', - packages=['multiprofit'], - # ext_modules=[], - # zip_safe=False, -) diff --git a/tests/SConscript b/tests/SConscript new file mode 100644 index 0000000..5437cb4 --- /dev/null +++ b/tests/SConscript @@ -0,0 +1,3 @@ +# -*- python -*- +from lsst.sconsUtils import scripts +scripts.BasicSConscript.tests(pyList=[]) diff --git a/tests/test_componentconfig.py b/tests/test_componentconfig.py index 772c233..fe8d284 100644 --- a/tests/test_componentconfig.py +++ b/tests/test_componentconfig.py @@ -1,5 +1,5 @@ -from multiprofit.config import set_config_from_dict -from multiprofit.componentconfig import EllipticalComponentConfig +from lsst.multiprofit.config import set_config_from_dict +from lsst.multiprofit.componentconfig import EllipticalComponentConfig def test_EllipticalComponentConfig(): diff --git a/tests/test_fit_bootstrap_model.py b/tests/test_fit_bootstrap_model.py index 2c03e21..2e3c832 100644 --- a/tests/test_fit_bootstrap_model.py +++ b/tests/test_fit_bootstrap_model.py @@ -20,17 +20,17 @@ # along with this program. If not, see . import gauss2d.fit as g2f -from multiprofit.componentconfig import ( +from lsst.multiprofit.componentconfig import ( GaussianConfig, init_component, ParameterConfig, SersicConfig, SersicIndexConfig, ) -from multiprofit.fit_bootstrap_model import ( +from lsst.multiprofit.fit_bootstrap_model import ( CatalogExposurePsfBootstrap, CatalogExposureSourcesBootstrap, CatalogSourceFitterBootstrap, ) -from multiprofit.fit_psf import CatalogPsfFitter, CatalogPsfFitterConfig -from multiprofit.fit_source import CatalogSourceFitterConfig -from multiprofit.modeller import ModelFitConfig -from multiprofit.plots import ErrorValues, plot_catalog_bootstrap, plot_loglike -from multiprofit.utils import get_params_uniq +from lsst.multiprofit.fit_psf import CatalogPsfFitter, CatalogPsfFitterConfig +from lsst.multiprofit.fit_source import CatalogSourceFitterConfig +from lsst.multiprofit.modeller import ModelFitConfig +from lsst.multiprofit.plots import ErrorValues, plot_catalog_bootstrap, plot_loglike +from lsst.multiprofit.utils import get_params_uniq import numpy as np import pytest diff --git a/tests/test_modeller.py b/tests/test_modeller.py index 17cc649..6ed642e 100644 --- a/tests/test_modeller.py +++ b/tests/test_modeller.py @@ -23,11 +23,11 @@ import gauss2d as g2 import gauss2d.fit as g2f import math -from multiprofit.modeller import ( +from lsst.multiprofit.modeller import ( fitmethods_linear, LinearGaussians, make_image_gaussians, make_psfmodel_null, Modeller, ) -from multiprofit.transforms import transforms_ref -from multiprofit.utils import get_params_uniq +from lsst.multiprofit.transforms import transforms_ref +from lsst.multiprofit.utils import get_params_uniq import numpy as np import pytest import scipy.optimize as spopt diff --git a/tests/test_psfmodel_utils.py b/tests/test_psfmodel_utils.py index db85936..1b09090 100644 --- a/tests/test_psfmodel_utils.py +++ b/tests/test_psfmodel_utils.py @@ -20,7 +20,7 @@ # along with this program. If not, see . import gauss2d.fit as g2f -from multiprofit.psfmodel_utils import make_psf_source +from lsst.multiprofit.psfmodel_utils import make_psf_source def test_make_psf_source(): diff --git a/ups/multiprofit.table b/ups/multiprofit.table new file mode 100644 index 0000000..0acbd3a --- /dev/null +++ b/ups/multiprofit.table @@ -0,0 +1,8 @@ +# List EUPS dependencies of this package here. +# - Any package whose API is used directly should be listed explicitly. +# - Common third-party packages can be assumed to be recursively included by +# the "sconsUtils" package. +setupRequired(sconsUtils) +setupRequired(gauss2dfit) + +envPrepend(PYTHONPATH, ${PRODUCT_DIR}/python)