From 2194f225c265902199db6e26a14d4b1fb48e6650 Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 14:29:04 -0600 Subject: [PATCH 1/6] Add fast build ellipse model --- photutils/isophote/model.py | 118 +++++++++---------------- photutils/isophote/tests/test_model.py | 88 ++++++++++++------ 2 files changed, 106 insertions(+), 100 deletions(-) diff --git a/photutils/isophote/model.py b/photutils/isophote/model.py index 407ae65a7..69bbe3b28 100644 --- a/photutils/isophote/model.py +++ b/photutils/isophote/model.py @@ -6,16 +6,15 @@ import numpy as np - -from .geometry import EllipseGeometry +import ctypes as ct +import os __all__ = ['build_ellipse_model'] -def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False): +def build_ellipse_model(shape, isolist, nthreads=1, fill=0., high_harmonics=False): """ Build a model elliptical galaxy image from a list of isophotes. - For each ellipse in the input isophote list the algorithm fills the output image array with the corresponding isophotal intensity. Pixels in the output array are in general only partially covered by @@ -24,27 +23,24 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False): each pixel by storing the partial area information in an auxiliary array. The information in this array is then used to normalize the pixel intensities. - Parameters ---------- shape : 2-tuple The (ny, nx) shape of the array used to generate the input ``isolist``. - isolist : `~photutils.isophote.IsophoteList` instance The isophote list created by the `~photutils.isophote.Ellipse` class. - + nthreads: float, optiomal + Number of threads to perform work. Default is 1 (serial code). fill : float, optional The constant value to fill empty pixels. If an output pixel has no contribution from any isophote, it will be assigned this value. The default is 0. - high_harmonics : bool, optional Whether to add the higher-order harmonics (i.e., ``a3``, ``b3``, ``a4``, and ``b4``; see `~photutils.isophote.Isophote` for details) to the result. - Returns ------- result : 2D `~numpy.ndarray` @@ -91,71 +87,42 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False): # correct deviations cased by fluctuations in spline solution eps_array[np.where(eps_array < 0.)] = 0. - - result = np.zeros(shape=shape) - weight = np.zeros(shape=shape) - eps_array[np.where(eps_array < 0.)] = 0.05 - - # for each interpolated isophote, generate intensity values on the - # output image array - # for index in range(len(finely_spaced_sma)): - for index in range(1, len(finely_spaced_sma)): - sma0 = finely_spaced_sma[index] - eps = eps_array[index] - pa = pa_array[index] - x0 = x0_array[index] - y0 = y0_array[index] - geometry = EllipseGeometry(x0, y0, sma0, eps, pa) - - intens = intens_array[index] - - # scan angles. Need to go a bit beyond full circle to ensure - # full coverage. - r = sma0 - phi = 0. - while phi <= 2*np.pi + geometry._phi_min: - # we might want to add the third and fourth harmonics - # to the basic isophotal intensity. - harm = 0. - if high_harmonics: - harm = (a3_array[index] * np.sin(3.*phi) + - b3_array[index] * np.cos(3.*phi) + - a4_array[index] * np.sin(4.*phi) + - b4_array[index] * np.cos(4.*phi)) / 4. - - # get image coordinates of (r, phi) pixel - x = r * np.cos(phi + pa) + x0 - y = r * np.sin(phi + pa) + y0 - i = int(x) - j = int(y) - - if (i > 0 and i < shape[1] - 1 and j > 0 and j < shape[0] - 1): - # get fractional deviations relative to target array - fx = x - float(i) - fy = y - float(j) - - # add up the isophote contribution to the overlapping pixels - result[j, i] += (intens + harm) * (1. - fy) * (1. - fx) - result[j, i + 1] += (intens + harm) * (1. - fy) * fx - result[j + 1, i] += (intens + harm) * fy * (1. - fx) - result[j + 1, i + 1] += (intens + harm) * fy * fx - - # add up the fractional area contribution to the - # overlapping pixels - weight[j, i] += (1. - fy) * (1. - fx) - weight[j, i + 1] += (1. - fy) * fx - weight[j + 1, i] += fy * (1. - fx) - weight[j + 1, i + 1] += fy * fx - - # step towards next pixel on ellipse - phi = max((phi + 0.75 / r), geometry._phi_min) - r = max(geometry.radius(phi), 0.5) - # if outside image boundaries, ignore. - else: - break - - # zero weight values must be set to 1. + + # convert everything to C-type array (pointers) + c_fss_array = ct.c_void_p(finely_spaced_sma.ctypes.data) + c_intens_array = ct.c_void_p(intens_array.ctypes.data) + c_eps_array = ct.c_void_p(eps_array.ctypes.data) + c_pa_array = ct.c_void_p(pa_array.ctypes.data) + c_x0_array = ct.c_void_p(x0_array.ctypes.data) + c_y0_array = ct.c_void_p(y0_array.ctypes.data) + c_a3_array = ct.c_void_p(a3_array.ctypes.data) + c_b3_array = ct.c_void_p(b3_array.ctypes.data) + c_a4_array = ct.c_void_p(a4_array.ctypes.data) + c_b4_array = ct.c_void_p(b4_array.ctypes.data) + + # initialize result and weight arrays, also as 1D ctype array + result = np.zeros(shape=(shape[1]*shape[0],)) + weight = np.zeros(shape=(shape[1]*shape[0],)) + c_result = ct.c_void_p(result.ctypes.data) + c_weight = ct.c_void_p(weight.ctypes.data) + + # convert high_harmnics bool flag to int, + # convert all other ints to ctype + c_high_harm = ct.c_int(int(high_harmonics)) + c_N = ct.c_int(len(finely_spaced_sma)) + c_nrows = ct.c_int(shape[0]) + c_ncols = ct.c_int(shape[1]) + + # load into C worker function (worker.so should be in same directory) + lib = ct.cdll.LoadLibrary(os.path.dirname(os.path.abspath(__file__)) + '/worker.so') + lib.worker.restype = None + lib.worker(c_result, c_weight, c_nrows, c_ncols, c_N, c_high_harm, + c_fss_array, c_intens_array, c_eps_array, c_pa_array, + c_x0_array, c_y0_array, c_a3_array, c_b3_array, + c_a4_array, c_b4_array, ct.c_int(nthreads)) + + # zero weight values must be set to 1. weight[np.where(weight <= 0.)] = 1. # normalize @@ -164,4 +131,7 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False): # fill value result[np.where(result == 0.)] = fill - return result + # reshape + result = result.reshape(shape[0], shape[1]) + + return result diff --git a/photutils/isophote/tests/test_model.py b/photutils/isophote/tests/test_model.py index 6484a992d..91d28b24d 100644 --- a/photutils/isophote/tests/test_model.py +++ b/photutils/isophote/tests/test_model.py @@ -2,12 +2,12 @@ """ Tests for the model module. """ -import warnings from astropy.io import fits -from astropy.utils.data import get_pkg_data_filename import numpy as np +import os.path as op import pytest +import multiprocessing as mp from .make_test_data import make_test_image from ..ellipse import Ellipse @@ -28,13 +28,7 @@ def test_model(): g = EllipseGeometry(530., 511, 10., 0.1, 10./180.*np.pi) ellipse = Ellipse(data, geometry=g, threshold=1.e5) - - # NOTE: this sometimes emits warnings (e.g., py38, ubuntu), but - # sometimes not. Here we simply ignore any RuntimeWarning, whether - # there is one or not. - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) - isophote_list = ellipse.fit_image() + isophote_list = ellipse.fit_image() model = build_ellipse_model(data.shape, isophote_list, fill=np.mean(data[10:100, 10:100])) @@ -43,7 +37,28 @@ def test_model(): residual = data - model assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 + +@pytest.mark.remote_data +@pytest.mark.skipif('not HAS_SCIPY') +def test_model_parallel(): + path = get_path('isophote/M105-S001-RGB.fits', + location='photutils-datasets', cache=True) + hdu = fits.open(path) + data = hdu[0].data[0] + hdu.close() + g = EllipseGeometry(530., 511, 10., 0.1, 10./180.*np.pi) + ellipse = Ellipse(data, geometry=g, threshold=1.e5) + isophote_list = ellipse.fit_image() + # test parallel on a thread pool, w/ size equal to CPU thread count + model = build_ellipse_model(data.shape, isophote_list, nthreads=mp.cpu_count(), + fill=np.mean(data[10:100, 10:100])) + + assert data.shape == model.shape + + residual = data - model + assert np.mean(residual) <= 5.0 + assert np.mean(residual) >= -5.0 @pytest.mark.skipif('not HAS_SCIPY') def test_model_simulated_data(): @@ -59,6 +74,26 @@ def test_model_simulated_data(): assert data.shape == model.shape residual = data - model + + assert np.mean(residual) <= 5.0 + assert np.mean(residual) >= -5.0 + +@pytest.mark.skipif('not HAS_SCIPY') +def test_model_large_array(): + data = make_test_image(nx=10000, ny=10000, i0=100., sma=500., eps=0.3, + pa=np.pi/4., noise=0.05, seed=0) + + g = EllipseGeometry(5000., 5000., 100., 0.3, np.pi/4.) + ellipse = Ellipse(data, geometry=g, threshold=1.e5) + isophote_list = ellipse.fit_image() + # test parallel on a thread pool, w/ size equal to CPU thread count + model = build_ellipse_model(data.shape, isophote_list, nthreads=mp.cpu_count(), + fill=np.mean(data[0:50, 0:50])) + + assert data.shape == model.shape + + residual = data - model + assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 @@ -69,20 +104,21 @@ def test_model_minimum_radius(): # model building algorithm into a corner, where it fails. # With the algorithm fixed, it bypasses the failure and # succeeds in building the model image. - filepath = get_pkg_data_filename('data/minimum_radius_test.fits') - with fits.open(filepath) as hdu: - data = hdu[0].data - - g = EllipseGeometry(50., 45, 530., 0.1, 10. / 180. * np.pi) - g.find_center(data) - ellipse = Ellipse(data, geometry=g) - with pytest.warns(RuntimeWarning, match='Degrees of freedom'): - isophote_list = ellipse.fit_image(sma0=40, minsma=0, maxsma=350., - step=0.4, nclip=3) - - model = build_ellipse_model(data.shape, isophote_list, - fill=np.mean(data[0:50, 0:50])) - - # It's enough that the algorithm reached this point. The - # actual accuracy of the modelling is being tested elsewhere. - assert data.shape == model.shape + filepath = op.join(op.dirname(op.abspath(__file__)), 'data', + 'minimum_radius_test.fits') + hdu = fits.open(filepath) + data = hdu[0].data + hdu.close() + + g = EllipseGeometry(50., 45, 530., 0.1, 10. / 180. * np.pi) + g.find_center(data) + ellipse = Ellipse(data, geometry=g) + isophote_list = ellipse.fit_image(sma0=40, minsma=0, maxsma=350., + step=0.4, nclip=3) + + model = build_ellipse_model(data.shape, isophote_list, + fill=np.mean(data[0:50, 0:50])) + + # It's enough that the algorithm reached this point. The + # actual accuracy of the modelling is being tested elsewhere. + assert data.shape == model.shape From b097aaab2c2ef6e8fb336640041803b47d0b985d Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 14:42:19 -0600 Subject: [PATCH 2/6] add shared object file --- photutils/isophote/worker.so | Bin 0 -> 17160 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100755 photutils/isophote/worker.so diff --git a/photutils/isophote/worker.so b/photutils/isophote/worker.so new file mode 100755 index 0000000000000000000000000000000000000000..290eea498d9bfacc242de13fc2feab5cf153cb34 GIT binary patch literal 17160 zcmeHOeQ+DcbzcyqM8}o@*>X$Ui3KZ?OS=|K`BNWNe6)DdIgxBxA04uxeS!#3WJMA* z09s<|y5um(5ChYtlXymU#xqJ=O_U}XE0a1^Khz!|_mOGE@JB<>20XzkToR+jo0=x5wK*GI#B02m}O|p!k>|X?~SL3|%p~ zRWl%l*d(sR@ka3xsoPkoTEa~?D~PEGE1BA8B38>@`b`4yD~y@z2Z^Qr6KkML>UtJ|w}9t7uKQ=ML! zPx3)yX6J(kNNt8<3F?N4U#(@TPYn&d_TAAZAN}26+gHzzum8z`@!%+_lH05CBfmQS zF|ZN*ByRuv)OzA;l{ZmA_TECVMZ6&7zroK>f3ODr$2IUD)xb~Gz`s-jr+!E9^Oi>d z@HukV9jJkK)xcvl@Y{e#@bi{&0Dk%TZVh|~@U>nikb9;6N-=zms{Cd^E7|tp!Q_yg z&csp~+ZOf%`*+*%L@IH(Kb=XW_V3<4kQ_?vj~yCFs5<`|yXRQU?(H9n4fKC1A$pQ& zF`PUq()~js{a7j^jwVw_5-E`!9JbSmj6F0mXlMFTiC8=>c0ROwk3Aeq#Rdiv1Hz_2 zdXCsVeMjuxSpR_N0iNuUmWMr-11>x2|dK!qqOgewkc*KF}ohPf4hB0^e%_HBYHQIM^{P^e~q~28l3A)Wp$avD8&KP zybVfY$G~syqc`JHTqeS(By6sgQmW_2J2gUXJ`~={8nXZGnSl!>hf@b)&8XLL3 z{|bubHm%R|HP4sRdhj&sd71G&8Xq&rYh`|AdHz0y(;PrO2-#!MQMLOa#mZx+p#;*uQnB)Q&`<*5S1DE=`wS%zjuYsVVMDQm>sz`6j=6BZ zc}H2*kM#`HGf>Y!Jp=U&)HCq^k%2ctH@;&{TyD1VP2ay=2rE02X(-QG6W?w=Cz~tZ z{UgAo^7{V-39S*N-$ym2;(wIOC7Dt_IO=>exy{Gub% zgH5UUj3d(nO{qBU$n@|~Dn91O^Z-*T#vGX*)JnyN9hn|vO2u|ZrU#f(aib&CgG;G+ z6U*q`PqDlYt?3Z5e~|yjBfsU5z51_t>i^9n|BFZdCy)HHM?UM3U-HOb^T;PX@(V0m zQ`^1>-kNH9p2+6K%(XbxKi;5D^%dp$(3&x1L!8qT`j1KdyMLd^E$>;m1#9+a_gS-x zK`U_HdgZ;$YB+dFJ7_M?_sYDnJ&lL4Z5v=KMt=K%HL>lp&=gkgwaiske%ngmh1`4P zav_f7pKoe{9C!q_PWxnEJc>$U4#2JyT4gTcuzR5t3Y&|)q1(-c$Cpd~Crv~a9y!Ol z7{~_c7f62|hXa?O(bR$m(o(3~gwSOYY|arJ0NfjD1<;*9)vClrD{mHhtF)Sm`lm&ViNfqVw&w!PY=FHPzW&Anw6j<4-oOuGQ z>}$?E4%R~BXkt2_-9y2cm-8}BIkS*g;hC?haC7E7)DJ54Nw9s29kZs)vA;s!#9!0F z%ALoM|1h}1*U_74BziG!-Kn~tX)1geVCAf}@84&QUpCZD*Kh8Wytej#Uc~ znF0q10~xMZ<&A+{wovIq-ssC^FT3SX)fKWKgCgKr;j;jaBAx9~GI$ipmKe6`A~s$|)ZedB!R#r+rl9X|AZ8@llauqoOiH zGuyrBH^ok~mw6nVyasziucLOAya=aJvzw}$u93Q(vzS6N9wvXTpPn@NZB0||w+{@ucnwCPH=IO$gx8>9+nI^RaH0 z^%R6-T~Q!g@R5z{H0Q%O2 z&}(?>r2gX?qnLm7#=e8Gz>R%N-coCN;0CDEwHrWB_-q@W2gF z(%!!T3YJo2&gW&CeJt4Fl8jAP^51j;gI)#Ot%Cigyb8K$DEX#a8ME&t>2~`*qFr&E zHI&)clDxQmt#O4#H^>^w>}yF9UA}7P&D}RTaet4eeLVT_^{=mb%r=Ev{<~bpT7bj@)Na-8|)$ z(uT49(hTb1V@fyIv`q%rHCy=fn^e)2RSK&l^>Md}#+7up3GU9G@wodLRyy!9TQIdd zQa}$d4mwY2*OL0UZgX6>!1W7W*I!`QKX^k$suVikFXC3mtxSN|Z3nykoUhyPijPW# zM`>n|7iIBA++hnRp?x7ckKU=fGI8pz{KVQqGb-`iG)edNzb#3An)uV;i`U>`DK9PE z=gAz?my=JGp6O~^@pz>m+)mZAhngmR{OZx{^nVH{8z zv!Iti+b~`hL4S;krUgl$LjSv^68O|!5jYkO+;G+M=4sTe14d8h$KZP&RICi|SlRON z(3MA<$HaZtZ@Kl(b!!PGdy1nCZMy%TKzKKPt-$9%xscBce$V~9T)vG!b7lB54ck{N z>p{a1ss8i$okaatK%M%AKamPl{}=eRp-=LEjcgzMT7VQHsQYL9&ZDo7()$Rm|M*kE ziN@e3KBSsKMZMHBP|rX;1N98lGf>Y!Jp=U&)H6`eKs^H=Oa}P-IsRTwFURCJeG~!9 z^o<@B{Hse9y~x7<&}8}a8bznJK*|JHluB*_bS4M7EO z0cDxTt1K$~ty>HIPZ@=lR#Z0W`1o7376n(nu`?8)-8%ZKwH~cssr;rQqEiEmpVkd? zoTs!NeiL@R#G%}PKiENJE8ExCJ@!kfsE-`PsO9hN-PJ*S>B-;ISK~ovyjA0kg2!1k-XwTDR^!VApVw-9x#00o zjW-J(SJn6x&iJdwG1_(Bs_`pD#23fa>i#HI8wK%&X4n^B>FoQc_8$^_eOKdXG~z}S zA4H>uFMhR{^u_VfHt)x9xeMa04$l)V7}xZ*!f)IL#cv3HLshMhtJ-iQiVunw!Q;%X#L_$;2ksaDmuldrfs_B|bo}z&88CVmM)=Dbm+#Dg(f$X* zZ@v}`ihLIaS}Sk=DC@2kPf$W!rBkSz&LQvP@@3F=YU5fSABS~ zhWFRrebFX?H<~wlEw~_c+x&RkUWG9An{Bxm9}Fe z$FL)1aCjh*NyMX@+U{xdSEBtZ{dO#sijCTdp-gHN8%1J+2|GSAI5-Lumt;dTQ>B{h z?Hx#BKZ_1cZFoV9gRx^G6^r+er1d74s7?D;dWW#H#I|?rZQpI$=7Sxy^M;xSrajq@ z#}GShSE)d~BSSrm#nVZP zvXj1hpQ<)6=}4RnL{5})f0ZjvyK2W&b-W^)9v#Gq0L`Qn?PG2znMp(s4~;~JQ_0~( zDl;mg4s>Xwe;~f0KQ5#~Uo71xqVdroI8!u}QZdi^CgB%7?mGlm+u>DmK^;XK+oP z5vBBVD4iBjT(N_=3jL;lUk9Ug7Ohn>{6Go}B-xr<(o zIZ!(|{J);2Q(1WZ#+2>Z&fEScU`1;;w&&+=rWbLap+f6HV%eV8UHgF18j9`t`J5@Q z51}H9cE_UlBsf~Du{}S>Grd`Ba{R2vv=8;P_G6i!_n8iBdy1cQIDQR_qhQE4+w;1B z=_#$q?Xx}ie_Y$&ul0ES%9Q(0_ITQG+GqPbFv^KumR4?}`ew@Oc&Nz2aj^UuFm8KZ zk1(BQMNRo0&GZ?My|Zr7>lD^wH!L%K-eb?}8>Vw?s43eq{Zo(qyl$82pYe+hk3GBn zvd5mEJDGaS<+)~CqAzNDKK~*@C`v>gg`tRZDq z8H{Or)>G)xe9{7xbU1!qXH3?xr?%Xs5ddLY-Kk>xTOC9(Z11K#_S~vxw8Kk!Vxjka p_u|L>;rZtZU>v*Kewvg~UiCQVc2%u<&EaRC*`hR7cnmzE_+Q!w=End4 literal 0 HcmV?d00001 From d6c310f71a687a90136f5be212be4f111d5e8e74 Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 15:00:25 -0600 Subject: [PATCH 3/6] add new tests for parallel and large-array --- photutils/isophote/tests/test_model.py | 58 ++++++++++++++------------ 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/photutils/isophote/tests/test_model.py b/photutils/isophote/tests/test_model.py index 91d28b24d..872ee1247 100644 --- a/photutils/isophote/tests/test_model.py +++ b/photutils/isophote/tests/test_model.py @@ -2,12 +2,12 @@ """ Tests for the model module. """ +import warnings from astropy.io import fits +from astropy.utils.data import get_pkg_data_filename import numpy as np -import os.path as op import pytest -import multiprocessing as mp from .make_test_data import make_test_image from ..ellipse import Ellipse @@ -28,7 +28,13 @@ def test_model(): g = EllipseGeometry(530., 511, 10., 0.1, 10./180.*np.pi) ellipse = Ellipse(data, geometry=g, threshold=1.e5) - isophote_list = ellipse.fit_image() + + # NOTE: this sometimes emits warnings (e.g., py38, ubuntu), but + # sometimes not. Here we simply ignore any RuntimeWarning, whether + # there is one or not. + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + isophote_list = ellipse.fit_image() model = build_ellipse_model(data.shape, isophote_list, fill=np.mean(data[10:100, 10:100])) @@ -37,7 +43,7 @@ def test_model(): residual = data - model assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 - + @pytest.mark.remote_data @pytest.mark.skipif('not HAS_SCIPY') def test_model_parallel(): @@ -60,6 +66,7 @@ def test_model_parallel(): assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 + @pytest.mark.skipif('not HAS_SCIPY') def test_model_simulated_data(): data = make_test_image(nx=200, ny=200, i0=10., sma=5., eps=0.5, @@ -74,16 +81,15 @@ def test_model_simulated_data(): assert data.shape == model.shape residual = data - model - assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 - + @pytest.mark.skipif('not HAS_SCIPY') def test_model_large_array(): - data = make_test_image(nx=10000, ny=10000, i0=100., sma=500., eps=0.3, + data = make_test_image(nx=5000, ny=5000, i0=100., sma=500., eps=0.3, pa=np.pi/4., noise=0.05, seed=0) - g = EllipseGeometry(5000., 5000., 100., 0.3, np.pi/4.) + g = EllipseGeometry(2500., 2500., 100., 0.3, np.pi/4.) ellipse = Ellipse(data, geometry=g, threshold=1.e5) isophote_list = ellipse.fit_image() # test parallel on a thread pool, w/ size equal to CPU thread count @@ -97,28 +103,26 @@ def test_model_large_array(): assert np.mean(residual) <= 5.0 assert np.mean(residual) >= -5.0 - @pytest.mark.skipif('not HAS_SCIPY') def test_model_minimum_radius(): # This test requires a "defective" image that drives the # model building algorithm into a corner, where it fails. # With the algorithm fixed, it bypasses the failure and # succeeds in building the model image. - filepath = op.join(op.dirname(op.abspath(__file__)), 'data', - 'minimum_radius_test.fits') - hdu = fits.open(filepath) - data = hdu[0].data - hdu.close() - - g = EllipseGeometry(50., 45, 530., 0.1, 10. / 180. * np.pi) - g.find_center(data) - ellipse = Ellipse(data, geometry=g) - isophote_list = ellipse.fit_image(sma0=40, minsma=0, maxsma=350., - step=0.4, nclip=3) - - model = build_ellipse_model(data.shape, isophote_list, - fill=np.mean(data[0:50, 0:50])) - - # It's enough that the algorithm reached this point. The - # actual accuracy of the modelling is being tested elsewhere. - assert data.shape == model.shape + filepath = get_pkg_data_filename('data/minimum_radius_test.fits') + with fits.open(filepath) as hdu: + data = hdu[0].data + + g = EllipseGeometry(50., 45, 530., 0.1, 10. / 180. * np.pi) + g.find_center(data) + ellipse = Ellipse(data, geometry=g) + with pytest.warns(RuntimeWarning, match='Degrees of freedom'): + isophote_list = ellipse.fit_image(sma0=40, minsma=0, maxsma=350., + step=0.4, nclip=3) + + model = build_ellipse_model(data.shape, isophote_list, + fill=np.mean(data[0:50, 0:50])) + + # It's enough that the algorithm reached this point. The + # actual accuracy of the modelling is being tested elsewhere. + assert data.shape == model.shape From ba94c602da2d62e60b47ad543dc42e9f3331835b Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 15:15:20 -0600 Subject: [PATCH 4/6] fix tests --- photutils/isophote/tests/test_model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/photutils/isophote/tests/test_model.py b/photutils/isophote/tests/test_model.py index 872ee1247..2e1bddba2 100644 --- a/photutils/isophote/tests/test_model.py +++ b/photutils/isophote/tests/test_model.py @@ -8,6 +8,7 @@ from astropy.utils.data import get_pkg_data_filename import numpy as np import pytest +import multiprocessing as mp from .make_test_data import make_test_image from ..ellipse import Ellipse From dc1d4a92a2bd2bee71d247da7e3007cee362aaa6 Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 15:24:14 -0600 Subject: [PATCH 5/6] fix warnings for tests --- photutils/isophote/tests/test_model.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/photutils/isophote/tests/test_model.py b/photutils/isophote/tests/test_model.py index 2e1bddba2..4df714520 100644 --- a/photutils/isophote/tests/test_model.py +++ b/photutils/isophote/tests/test_model.py @@ -56,7 +56,13 @@ def test_model_parallel(): g = EllipseGeometry(530., 511, 10., 0.1, 10./180.*np.pi) ellipse = Ellipse(data, geometry=g, threshold=1.e5) - isophote_list = ellipse.fit_image() + # NOTE: this sometimes emits warnings (e.g., py38, ubuntu), but + # sometimes not. Here we simply ignore any RuntimeWarning, whether + # there is one or not. + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + isophote_list = ellipse.fit_image() + # test parallel on a thread pool, w/ size equal to CPU thread count model = build_ellipse_model(data.shape, isophote_list, nthreads=mp.cpu_count(), fill=np.mean(data[10:100, 10:100])) From 97af14d2692cd0a541e00a48798af3648672e0eb Mon Sep 17 00:00:00 2001 From: Huan Le Date: Mon, 31 Jan 2022 15:37:45 -0600 Subject: [PATCH 6/6] binary file from fixed C code --- photutils/isophote/worker.so | Bin 17160 -> 17160 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/photutils/isophote/worker.so b/photutils/isophote/worker.so index 290eea498d9bfacc242de13fc2feab5cf153cb34..aeb34aa3d613ca49d25c420a74930d8f9b47327e 100755 GIT binary patch delta 530 zcmeBZW9(>SoG^n^m!AO)bT=+EXA<#xyDY=viq{NRHIuxE+NC=apRS18{Dx_zuz*YV z(l`J9|9^4**Z=>WhhJ1|J}4%_$QV5Nnz%Nj{$wVJ0LHM%!4k@hN}J0i{27IU*ccc- z@yEEn_{1N1;S;|AtK;MwlF^J$C!0v6GS*L?B{f_4F$)7j8ovfmHYOCPOn~(g%jBDq za+4F}Ste&mt8=yk<-tZwc94>QigHapAT7&iJ^7w=iS#pOkbGRL!2kdMKk-MjN`Pqr z-xDufC+EuO@jHD3i3{+~1*&_&`EBwl8CQWwu%5_+FBW_ONk&fID6248URIuq3nq~` z*+VvvYdaIjY@oVtlNZbSS|q@@>pz3FweA201<1|=U|N9J9;_hl;3s~;Rt2CtUdjQ3 z>BYgx8geoUDnLCR-L4%ToyT9?_ylH80I|;i*)LBrFif5$CpEcFZiBQ(;}Hi&28L+I znApP%44?Vset9&%;qYia!f|-=M2?%A*%aa!C4O)Mqsx^~ppD6ym(8n(xtWFeDJKs* z$L2=GdY;JvOstdiJQx)=d+Ob122z?v?o3QSI44^f%T1nO5&#r=XyVSrXg68YUY}Dz O9uoBmn-lGqGXnt1D$>FL delta 569 zcmeBZW9(>SoG^oPB0mEdOx(E8oJpi8FYoEO%DT8PW`i}y8dYD1HZoUke#5j8cd*<3=Ik!^Fdgg>K@85;w`C;k}M7oYeeFMQ$`U{#&GM>3l6wnTEO?j3(d)aGJ5=K-$3F5yp2G0FW!Eg+$ZD8Wee65dGJNcG721@ z|Nr;scJ1)!JpN+OCop>gh`j;GekshzFu6`nYI2g?hA>DVdd%?X6*&qDRInz-3u*lC z4*<=DhB1ot1P_FCGdKzuU#foxTG4pKfsuhB+A${f@MHz~rp;^$af}jiTtGuz`2^aS zoO#*2dYGG8n7z4p*f};gD%SH%4q#%Pr02odve{GbJ~NQgG;(KRisPDWWh^&&f=K{S d7vqx2n)dpfE%G1-GBC7kPPAXn3;>l$