From 9ecbad08dbecff2ac41fcd31639f027769167aa7 Mon Sep 17 00:00:00 2001 From: PaulHancock Date: Mon, 15 Nov 2021 14:06:16 +0800 Subject: [PATCH 1/3] pep8 and small tidyup --- AegeanTools/fitting.py | 164 ++++++++++++++++++++++++++--------------- tests/test_fitting.py | 64 +++++++++++----- 2 files changed, 148 insertions(+), 80 deletions(-) diff --git a/AegeanTools/fitting.py b/AegeanTools/fitting.py index 82b5d099..ac7e29e9 100644 --- a/AegeanTools/fitting.py +++ b/AegeanTools/fitting.py @@ -58,7 +58,7 @@ def elliptical_gaussian(x, y, amp, xo, yo, sx, sy, theta): xxo = x - xo yyo = y - yo exp = (xxo * cost + yyo * sint) ** 2 / sx ** 2 \ - + (xxo * sint - yyo * cost) ** 2 / sy ** 2 + + (xxo * sint - yyo * cost) ** 2 / sy ** 2 exp *= -1. / 2 return amp * np.exp(exp) @@ -100,7 +100,7 @@ def elliptical_gaussian_with_alpha(x, y, v, amp, xo, yo, vo, sx, sy, theta, alph if beta is not None: exponent += beta * np.log10(v/vo) snu = amp * (v/vo) ** (exponent) - gauss = elliptical_gaussian(x,y,snu,xo,yo,sx,sy,theta) + gauss = elliptical_gaussian(x, y, snu, xo, yo, sx, sy, theta) return gauss @@ -124,7 +124,8 @@ def Cmatrix(x, y, sx, sy, theta): data : array-like The C-matrix. """ - C = np.vstack([elliptical_gaussian(x, y, 1, i, j, sx, sy, theta) for i, j in zip(x, y)]) + C = np.vstack([elliptical_gaussian(x, y, 1, i, j, sx, sy, theta) + for i, j in zip(x, y)]) return C @@ -205,12 +206,14 @@ def jacobian(pars, x, y): matrix.append(dmds) if pars[prefix + 'xo'].vary: - dmdxo = cost * (xcos + ysin) / sx ** 2 + sint * (xsin - ycos) / sy ** 2 + dmdxo = cost * (xcos + ysin) / sx ** 2 + \ + sint * (xsin - ycos) / sy ** 2 dmdxo *= model matrix.append(dmdxo) if pars[prefix + 'yo'].vary: - dmdyo = sint * (xcos + ysin) / sx ** 2 - cost * (xsin - ycos) / sy ** 2 + dmdyo = sint * (xcos + ysin) / sx ** 2 - \ + cost * (xsin - ycos) / sy ** 2 dmdyo *= model matrix.append(dmdyo) @@ -223,7 +226,8 @@ def jacobian(pars, x, y): matrix.append(dmdsy) if pars[prefix + 'theta'].vary: - dmdtheta = model * (sy ** 2 - sx ** 2) * (xsin - ycos) * (xcos + ysin) / sx ** 2 / sy ** 2 + dmdtheta = model * (sy ** 2 - sx ** 2) * \ + (xsin - ycos) * (xcos + ysin) / sx ** 2 / sy ** 2 matrix.append(dmdtheta) return np.array(matrix) @@ -428,27 +432,30 @@ def hessian(pars, x, y): # H(xo,xo)/G = 1.0*(-sx**2*sy**2*(sx**2*sin(t)**2 + sy**2*cos(t)**2) + (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t))**2)/(sx**4*sy**4) hmat[j][k] = -sx2*sy2*(sx2*sint**2 + sy2*cost**2) hmat[j][k] += (sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost)**2 - hmat[j][k] *= model/ (sx2**2*sy2**2) + hmat[j][k] *= model / (sx2**2*sy2**2) k += 1 if yo_var: # H(xo,yo)/G = 1.0*(sx**2*sy**2*(sx**2 - sy**2)*sin(2*t)/2 - (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**4*sy**4) hmat[j][k] = sx2*sy2*(sx2 - sy2)*sin2t/2 - hmat[j][k] -= (sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost)*(sx2*(xsin -ycos)*cost - sy2*(xcos + ysin)*sint) + hmat[j][k] -= (sx2*(xsin - ycos)*sint + sy2*(xcos + ysin) + * cost)*(sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) hmat[j][k] *= model / (sx**4*sy**4) k += 1 if sx_var: # H(xo,sx) = ((x - xo)*cos(t) + (y - yo)*sin(t))*(-2.0*sx**2*sy**2*cos(t) + 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx**5*sy**2) hmat[j][k] = (xcos + ysin) - hmat[j][k] *= -2*sx2*sy2*cost + (xcos + ysin)*(sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost) + hmat[j][k] *= -2*sx2*sy2*cost + \ + (xcos + ysin)*(sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost) hmat[j][k] *= model / (sx**5*sy2) k += 1 if sy_var: # H(xo,sy) = ((x - xo)*sin(t) + (-y + yo)*cos(t))*(-2.0*sx**2*sy**2*sin(t) + 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx2*sy**5) hmat[j][k] = (xsin - ycos) - hmat[j][k] *= -2*sx2*sy2*sint + (xsin - ycos)*(sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost) + hmat[j][k] *= -2*sx2*sy2*sint + \ + (xsin - ycos)*(sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost) hmat[j][k] *= model/(sx2*sy**5) k += 1 @@ -456,9 +463,9 @@ def hessian(pars, x, y): # H(xo,t) = 1.0*(sx**2*sy**2*(sx**2 - sy**2)*(x*sin(2*t) - xo*sin(2*t) - y*cos(2*t) + yo*cos(2*t)) + (-sx**2 + 1.0*sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx**4*sy**4) # second part hmat[j][k] = (sy2-sx2)*(xsin - ycos)*(xcos + ysin) - hmat[j][k] *= sx2*(xsin -ycos)*sint + sy2*(xcos + ysin)*cost + hmat[j][k] *= sx2*(xsin - ycos)*sint + sy2*(xcos + ysin)*cost # first part - hmat[j][k] += sx2*sy2*(sx2 - sy2)*(xxo*sin2t -yyo*cos2t) + hmat[j][k] += sx2*sy2*(sx2 - sy2)*(xxo*sin2t - yyo*cos2t) hmat[j][k] *= model/(sx**4*sy**4) # k += 1 j += 1 @@ -472,12 +479,13 @@ def hessian(pars, x, y): if xo_var: # H(yo,xo)/G = H(xo,yo)/G - hmat[j][k] =hmat[1][2] + hmat[j][k] = hmat[1][2] k += 1 # if yo_var: # H(yo,yo)/G = 1.0*(-sx**2*sy**2*(sx**2*cos(t)**2 + sy**2*sin(t)**2) + (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t))**2)/(sx**4*sy**4) - hmat[j][k] = (sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint)**2 / (sx2**2*sy2**2) + hmat[j][k] = (sx2*(xsin - ycos)*cost - sy2 * + (xcos + ysin)*sint)**2 / (sx2**2*sy2**2) hmat[j][k] -= cost**2/sy2 + sint**2/sx2 hmat[j][k] *= model k += 1 @@ -485,14 +493,16 @@ def hessian(pars, x, y): if sx_var: # H(yo,sx)/G = -((x - xo)*cos(t) + (y - yo)*sin(t))*(2.0*sx**2*sy**2*sin(t) + 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) - (y - yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**5*sy**2) hmat[j][k] = -1*(xcos + ysin) - hmat[j][k] *= 2*sx2*sy2*sint + (xcos + ysin)*(sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) + hmat[j][k] *= 2*sx2*sy2*sint + \ + (xcos + ysin)*(sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) hmat[j][k] *= model/(sx**5*sy2) k += 1 if sy_var: # H(yo,sy)/G = ((x - xo)*sin(t) + (-y + yo)*cos(t))*(2.0*sx**2*sy**2*cos(t) - 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**2*sy**5) - hmat[j][k] = (xsin -ycos) - hmat[j][k] *= 2*sx2*sy2*cost - (xsin - ycos)*(sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) + hmat[j][k] = (xsin - ycos) + hmat[j][k] *= 2*sx2*sy2*cost - \ + (xsin - ycos)*(sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) hmat[j][k] *= model/(sx2*sy**5) k += 1 @@ -500,7 +510,8 @@ def hessian(pars, x, y): # H(yo,t)/G = 1.0*(sx**2*sy**2*(sx**2*(-x*cos(2*t) + xo*cos(2*t) - y*sin(2*t) + yo*sin(2*t)) + sy**2*(x*cos(2*t) - xo*cos(2*t) + y*sin(2*t) - yo*sin(2*t))) + (1.0*sx**2 - sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**4*sy**4) hmat[j][k] = (sx2 - sy2)*(xsin - ycos)*(xcos + ysin) hmat[j][k] *= (sx2*(xsin - ycos)*cost - sy2*(xcos + ysin)*sint) - hmat[j][k] += sx2*sy2*(sx2-sy2)*(-x*cos2t + xo*cos2t - y*sin2t + yo*sin2t) + hmat[j][k] += sx2*sy2 * \ + (sx2-sy2)*(-x*cos2t + xo*cos2t - y*sin2t + yo*sin2t) hmat[j][k] *= model/(sx**4*sy**4) # k += 1 j += 1 @@ -538,7 +549,7 @@ def hessian(pars, x, y): if theta_var: # H(sx,t)/G = (-2.0*sx**2*sy**2 + 1.0*(-sx**2 + sy**2)*((x - xo)*cos(t) + (y - yo)*sin(t))**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))/(sx**5*sy**2) hmat[j][k] = -2*sx2*sy2 + (sy2 - sx2)*(xcos + ysin)**2 - hmat[j][k] *= (xsin -ycos)*(xcos + ysin) + hmat[j][k] *= (xsin - ycos)*(xcos + ysin) hmat[j][k] *= model/(sx**5*sy**2) # k += 1 j += 1 @@ -602,7 +613,8 @@ def hessian(pars, x, y): # if theta_var: # H(t,t)/G = (sx**2*sy**2*(sx**2*(((x - xo)*sin(t) + (-y + yo)*cos(t))**2 - 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))**2) + sy**2*(-1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))**2 + ((x - xo)*cos(t) + (y - yo)*sin(t))**2)) + (sx**2 - 1.0*sy**2)**2*((x - xo)*sin(t) + (-y + yo)*cos(t))**2*((x - xo)*cos(t) + (y - yo)*sin(t))**2)/(sx**4*sy**4) hmat[j][k] = sx2*sy2 - hmat[j][k] *= sx2*((xsin - ycos)**2 - (xcos + ysin)**2) + sy2*((xcos + ysin)**2 - (xsin - ycos)**2) + hmat[j][k] *= sx2*((xsin - ycos)**2 - (xcos + ysin)**2) + \ + sy2*((xcos + ysin)**2 - (xsin - ycos)**2) hmat[j][k] += (sx2 - sy2)**2*(xsin - ycos)**2*(xcos + ysin)**2 hmat[j][k] *= model/(sx**4*sy**4) # j += 1 @@ -673,7 +685,7 @@ def nan_acf(noise): The ACF. """ corr = np.zeros(noise.shape) - ix,jx = noise.shape + ix, jx = noise.shape for i in range(ix): si_min = slice(i, None, None) si_max = slice(None, ix-i, None) @@ -683,7 +695,8 @@ def nan_acf(noise): if np.all(np.isnan(noise[si_min, sj_min])) or np.all(np.isnan(noise[si_max, sj_max])): corr[i, j] = np.nan else: - corr[i, j] = np.nansum(noise[si_min, sj_min] * noise[si_max, sj_max]) + corr[i, j] = np.nansum( + noise[si_min, sj_min] * noise[si_max, sj_max]) # return the normalised acf return corr / np.nanmax(corr) @@ -756,7 +769,8 @@ def RB_bias(data, pars, ita=None, acf=None): # all pixels x, y = np.indices(data.shape) # Create the jacobian as an AxN array accounting for the masked pixels - j = np.array(np.vsplit(lmfit_jacobian(pars, xm, ym).T, nparams)).reshape(nparams, -1) + j = np.array(np.vsplit(lmfit_jacobian(pars, xm, ym).T, + nparams)).reshape(nparams, -1) h = hessian(pars, x, y) # mask the hessian to be AxAxN array @@ -766,7 +780,7 @@ def RB_bias(data, pars, ita=None, acf=None): Bijk = np.einsum('ip,jkp', j, h) Eilkm = np.einsum('il,km', Dij, Dij) - Cimn_1 = -1 * np.einsum('krj,ir,km,jn', Bijk, Dij, Dij, Dij) + Cimn_1 = -1 * np.einsum('krj,ir,km,jn', Bijk, Dij, Dij, Dij) Cimn_2 = -1./2 * np.einsum('rkj,ir,km,jn', Bijk, Dij, Dij, Dij) Cimn = Cimn_1 + Cimn_2 @@ -878,11 +892,12 @@ def condon_errors(source, theta_n, psf=None): return smoothing = major * minor / (theta_n ** 2) - factor1 = (1 + (theta_n/ major)**2) - factor2 = (1 + (theta_n/ minor)**2) + factor1 = (1 + (theta_n / major)**2) + factor2 = (1 + (theta_n / minor)**2) snr = source.peak_flux / source.local_rms # calculation of rho2 depends on the parameter being used so we lambda this into a function - rho2 = lambda x: smoothing / 4 * factor1 ** alphas[x][0] * factor2 ** alphas[x][1] * snr ** 2 + def rho2(x): return smoothing / 4 * \ + factor1 ** alphas[x][0] * factor2 ** alphas[x][1] * snr ** 2 source.err_peak_flux = source.peak_flux * np.sqrt(2 / rho2('amp')) source.err_a = major * np.sqrt(2 / rho2('major')) * 3600. # arcsec @@ -890,8 +905,10 @@ def condon_errors(source, theta_n, psf=None): err_xo2 = 2. / rho2('xo') * major ** 2 / (8 * np.log(2)) # Condon'97 eq 21 err_yo2 = 2. / rho2('yo') * minor ** 2 / (8 * np.log(2)) - source.err_ra = np.sqrt(err_xo2 * np.sin(phi)**2 + err_yo2 * np.cos(phi)**2) - source.err_dec = np.sqrt(err_xo2 * np.cos(phi)**2 + err_yo2 * np.sin(phi)**2) + source.err_ra = np.sqrt(err_xo2 * np.sin(phi) ** + 2 + err_yo2 * np.cos(phi)**2) + source.err_dec = np.sqrt(err_xo2 * np.cos(phi) ** + 2 + err_yo2 * np.sin(phi)**2) if (major == 0) or (minor == 0): source.err_pa = ERR_MASK @@ -899,11 +916,13 @@ def condon_errors(source, theta_n, psf=None): elif abs(2 * (major-minor) / (major+minor)) < 0.01: source.err_pa = ERR_MASK else: - source.err_pa = np.degrees(np.sqrt(4 / rho2('pa')) * (major * minor / (major ** 2 - minor ** 2))) + source.err_pa = np.degrees( + np.sqrt(4 / rho2('pa')) * (major * minor / (major ** 2 - minor ** 2))) # integrated flux error err2 = (source.err_peak_flux / source.peak_flux) ** 2 - err2 += (theta_n ** 2 / (major * minor)) * ((source.err_a / source.a) ** 2 + (source.err_b / source.b) ** 2) + err2 += (theta_n ** 2 / (major * minor)) * \ + ((source.err_a / source.a) ** 2 + (source.err_b / source.b) ** 2) source.err_int_flux = source.int_flux * np.sqrt(err2) return @@ -974,23 +993,27 @@ def errors(source, model, wcshelper): if model[prefix + 'theta'].vary and np.isfinite(err_theta): # pa error - off1 = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) + off1 = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) off2 = wcshelper.pix2sky( [xo + sx * np.cos(np.radians(theta + err_theta)), yo + sy * np.sin(np.radians(theta + err_theta))]) - source.err_pa = abs(bear(ref[0], ref[1], off1[0], off1[1]) - bear(ref[0], ref[1], off2[0], off2[1])) + source.err_pa = abs( + bear(ref[0], ref[1], off1[0], off1[1]) - bear(ref[0], ref[1], off2[0], off2[1])) else: source.err_pa = ERR_MASK if model[prefix + 'sx'].vary and model[prefix + 'sy'].vary \ and all(np.isfinite([err_sx, err_sy])): # major axis error - ref = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) + ref = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) offset = wcshelper.pix2sky( [xo + (sx + err_sx) * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) source.err_a = gcd(ref[0], ref[1], offset[0], offset[1]) * 3600 # minor axis error - ref = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta + 90)), yo + sy * np.sin(np.radians(theta + 90))]) + ref = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta + 90)), yo + sy * np.sin(np.radians(theta + 90))]) offset = wcshelper.pix2sky( [xo + sx * np.cos(np.radians(theta + 90)), yo + (sy + err_sy) * np.sin(np.radians(theta + 90))]) source.err_b = gcd(ref[0], ref[1], offset[0], offset[1]) * 3600 @@ -998,7 +1021,8 @@ def errors(source, model, wcshelper): source.err_a = source.err_b = ERR_MASK sqerr = 0 - sqerr += (source.err_peak_flux / source.peak_flux) ** 2 if source.err_peak_flux > 0 else 0 + sqerr += (source.err_peak_flux / + source.peak_flux) ** 2 if source.err_peak_flux > 0 else 0 sqerr += (source.err_a / source.a) ** 2 if source.err_a > 0 else 0 sqerr += (source.err_b / source.b) ** 2 if source.err_b > 0 else 0 if sqerr == 0: @@ -1084,43 +1108,55 @@ def new_errors(source, model, wcshelper): # pragma: no cover (a, b), e = np.linalg.eig(mat) pa = np.degrees(np.arctan2(*e[0])) # transform this ellipse into sky coordinates - _, _, major, minor, pa = wcshelper.pix2sky_ellipse([xo, yo], a, b, pa) + _, _, major, minor, pa = wcshelper.pix2sky_ellipse( + [xo, yo], a, b, pa) # now determine the radius of the ellipse along the ra/dec directions. - source.err_ra = major*minor / np.hypot(major*np.sin(np.radians(pa)), minor*np.cos(np.radians(pa))) - source.err_dec = major*minor / np.hypot(major*np.cos(np.radians(pa)), minor*np.sin(np.radians(pa))) + source.err_ra = major*minor / \ + np.hypot(major*np.sin(np.radians(pa)), + minor*np.cos(np.radians(pa))) + source.err_dec = major*minor / \ + np.hypot(major*np.cos(np.radians(pa)), + minor*np.sin(np.radians(pa))) else: source.err_ra = source.err_dec = -1 if model[prefix + 'theta'].vary: # pa error - off1 = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) + off1 = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) # offset by 1 degree off2 = wcshelper.pix2sky( [xo + sx * np.cos(np.radians(theta + 1)), yo + sy * np.sin(np.radians(theta + 1))]) # scale the initial theta error by this amount - source.err_pa = abs(bear(ref[0], ref[1], off1[0], off1[1]) - bear(ref[0], ref[1], off2[0], off2[1])) * err_theta + source.err_pa = abs(bear(ref[0], ref[1], off1[0], off1[1]) - + bear(ref[0], ref[1], off2[0], off2[1])) * err_theta else: source.err_pa = ERR_MASK if model[prefix + 'sx'].vary and model[prefix + 'sy'].vary: # major axis error - ref = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) + ref = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) # offset by 0.1 pixels offset = wcshelper.pix2sky( [xo + (sx + 0.1) * np.cos(np.radians(theta)), yo + sy * np.sin(np.radians(theta))]) - source.err_a = gcd(ref[0], ref[1], offset[0], offset[1])/0.1 * err_sx * 3600 + source.err_a = gcd(ref[0], ref[1], offset[0], + offset[1])/0.1 * err_sx * 3600 # minor axis error - ref = wcshelper.pix2sky([xo + sx * np.cos(np.radians(theta + 90)), yo + sy * np.sin(np.radians(theta + 90))]) + ref = wcshelper.pix2sky( + [xo + sx * np.cos(np.radians(theta + 90)), yo + sy * np.sin(np.radians(theta + 90))]) # offset by 0.1 pixels offset = wcshelper.pix2sky( [xo + sx * np.cos(np.radians(theta + 90)), yo + (sy + 0.1) * np.sin(np.radians(theta + 90))]) - source.err_b = gcd(ref[0], ref[1], offset[0], offset[1])/0.1*err_sy * 3600 + source.err_b = gcd(ref[0], ref[1], offset[0], + offset[1])/0.1*err_sy * 3600 else: source.err_a = source.err_b = ERR_MASK sqerr = 0 - sqerr += (source.err_peak_flux / source.peak_flux) ** 2 if source.err_peak_flux > 0 else 0 + sqerr += (source.err_peak_flux / + source.peak_flux) ** 2 if source.err_peak_flux > 0 else 0 sqerr += (source.err_a / source.a) ** 2 if source.err_a > 0 else 0 sqerr += (source.err_b / source.b) ** 2 if source.err_b > 0 else 0 source.err_int_flux = abs(source.int_flux * np.sqrt(sqerr)) @@ -1240,9 +1276,11 @@ def residual(params, **kwargs): return (model - data[mask]).dot(B) if dojac: - result = lmfit.minimize(residual, params, kws={'x': mask[0], 'y': mask[1], 'B': B, 'errs': errs}, Dfun=lmfit_jacobian) + result = lmfit.minimize(residual, params, kws={ + 'x': mask[0], 'y': mask[1], 'B': B, 'errs': errs}, Dfun=lmfit_jacobian) else: - result = lmfit.minimize(residual, params, kws={'x': mask[0], 'y': mask[1], 'B': B, 'errs': errs}) + result = lmfit.minimize(residual, params, kws={ + 'x': mask[0], 'y': mask[1], 'B': B, 'errs': errs}) # Remake the residual so that it is once again (model - data) if B is not None: @@ -1311,7 +1349,6 @@ def covar_errors(params, data, errs, B, C=None): if __name__ == "__main__": - def plot_jacobian(): """ Plot the Jacobian for a test model @@ -1328,11 +1365,13 @@ def plot_jacobian(): # The model parameters params = lmfit.Parameters() params.add('c0_amp', value=1, min=0.5, max=2) - params.add('c0_xo', value=1. * nx / 2, min=nx / 2. - smoothing / 2., max=nx / 2. + smoothing / 2) - params.add('c0_yo', value=1. * ny / 2, min=ny / 2. - smoothing / 2., max=ny / 2. + smoothing / 2.) + params.add('c0_xo', value=1. * nx / 2, min=nx / 2. - + smoothing / 2., max=nx / 2. + smoothing / 2) + params.add('c0_yo', value=1. * ny / 2, min=ny / 2. - + smoothing / 2., max=ny / 2. + smoothing / 2.) params.add('c0_sx', value=2 * smoothing, min=0.8 * smoothing) params.add('c0_sy', value=smoothing, min=0.8 * smoothing) - params.add('c0_theta', value=45) #, min=-2*np.pi, max=2*np.pi) + params.add('c0_theta', value=45) # , min=-2*np.pi, max=2*np.pi) params.add('components', value=1, vary=False) def rmlabels(ax): @@ -1372,7 +1411,6 @@ def clx(ax): ax.set_yticks([]) return - def test_hessian_plots(): """ Plot the empirical and analytical hessian to check for agreement. @@ -1389,8 +1427,10 @@ def test_hessian_plots(): model.add('components', 1, vary=False) x, y = np.indices((40, 40)) # Empirical Hessian - kwargs = {"interpolation": "nearest", 'aspect': 1, 'vmin': -1, 'vmax': 1} - fig, ax = pyplot.subplots(6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) + kwargs = {"interpolation": "nearest", + 'aspect': 1, 'vmin': -1, 'vmax': 1} + fig, ax = pyplot.subplots( + 6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) Hemp = emp_hessian(model, x, y) params = ['amp', 'xo', 'yo', 'sx', 'sy', 'theta'] for i, row in enumerate(ax): @@ -1408,7 +1448,8 @@ def test_hessian_plots(): fig.suptitle('Empirical Hessian') # Analytical Hessian - fig, ax = pyplot.subplots(6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) + fig, ax = pyplot.subplots( + 6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) Hana = hessian(model, x, y) for i, row in enumerate(ax): for j, ax in enumerate(row): @@ -1425,7 +1466,8 @@ def test_hessian_plots(): fig.suptitle('Analytical Hessian') # Difference - fig, ax = pyplot.subplots(6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) + fig, ax = pyplot.subplots( + 6, 6, squeeze=True, sharex=True, sharey=True, figsize=(5, 6)) Hana = hessian(model, x, y) for i, row in enumerate(ax): for j, ax in enumerate(row): @@ -1442,7 +1484,6 @@ def test_hessian_plots(): fig.suptitle('Difference') pyplot.show() - def test_jacobian_plot(): """ @@ -1459,10 +1500,12 @@ def test_jacobian_plot(): model.add('components', 1, vary=False) x, y = np.indices((40, 40)) - kwargs = {"interpolation": "nearest", 'aspect': 1, 'vmin': -1, 'vmax': 1} + kwargs = {"interpolation": "nearest", + 'aspect': 1, 'vmin': -1, 'vmax': 1} var_names = ['amp', 'xo', 'yo', 'sx', 'sy', 'theta'] - fig, ax = pyplot.subplots(6, 3, sharex=True, sharey=True, figsize=(3, 6)) + fig, ax = pyplot.subplots( + 6, 3, sharex=True, sharey=True, figsize=(3, 6)) Jemp = emp_jacobian(model, x, y) Jana = jacobian(model, x, y) @@ -1487,4 +1530,3 @@ def test_jacobian_plot(): test_hessian_plots() test_jacobian_plot() - diff --git a/tests/test_fitting.py b/tests/test_fitting.py index 917f9a83..5f71aca5 100755 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -28,19 +28,25 @@ def make_model(): def test_elliptical_gaussian(): """Test our elliptical gaussian creation function""" x, y = np.indices((3, 3)) - gauss = fitting.elliptical_gaussian(x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=0) - if np.any(np.isnan(gauss)): raise AssertionError() - gauss = fitting.elliptical_gaussian(x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=np.inf) - if not (np.all(np.isnan(gauss))): raise AssertionError() + gauss = fitting.elliptical_gaussian( + x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=0) + if np.any(np.isnan(gauss)): + raise AssertionError() + gauss = fitting.elliptical_gaussian( + x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=np.inf) + if not (np.all(np.isnan(gauss))): + raise AssertionError() def test_CandBmatrix(): """Test that C and B matricies can be created without error""" x, y = map(np.ravel, np.indices((3, 3))) C = fitting.Cmatrix(x, y, sx=1, sy=2, theta=0) - if np.any(np.isnan(C)): raise AssertionError() + if np.any(np.isnan(C)): + raise AssertionError() B = fitting.Bmatrix(C) - if np.any(np.isnan(B)): raise AssertionError() + if np.any(np.isnan(B)): + raise AssertionError() def test_hessian_shape(): @@ -50,7 +56,8 @@ def test_hessian_shape(): nvar = 3 x, y = np.indices((10, 10)) Hij = fitting.hessian(model, x, y) - if not (Hij.shape == (nvar, nvar, 10, 10)): raise AssertionError() + if not (Hij.shape == (nvar, nvar, 10, 10)): + raise AssertionError() # now add another component model.add('c1_amp', 1, vary=True) @@ -62,7 +69,8 @@ def test_hessian_shape(): nvar = 9 model['components'].value = 2 Hij = fitting.hessian(model, x, y) - if not (Hij.shape == (nvar, nvar, 10, 10)): raise AssertionError() + if not (Hij.shape == (nvar, nvar, 10, 10)): + raise AssertionError() def test_jacobian_shape(): @@ -74,7 +82,8 @@ def test_jacobian_shape(): nvar = 3 x, y = np.indices((10, 10)) Jij = fitting.jacobian(model, x, y) - if not (Jij.shape == (nvar, 10, 10)): raise AssertionError() + if not (Jij.shape == (nvar, 10, 10)): + raise AssertionError() model.add('c1_amp', 1, vary=True) model.add('c1_xo', 5, vary=True) @@ -85,7 +94,8 @@ def test_jacobian_shape(): nvar = 9 model['components'].value = 2 Jij = fitting.jacobian(model, x, y) - if not (Jij.shape == (nvar, 10, 10)): raise AssertionError() + if not (Jij.shape == (nvar, 10, 10)): + raise AssertionError() def test_emp_vs_ana_jacobian(): @@ -96,7 +106,8 @@ def test_emp_vs_ana_jacobian(): emp_Jij = fitting.emp_jacobian(model, x, y) ana_Jij = fitting.jacobian(model, x, y) diff = np.abs(ana_Jij - emp_Jij) - if not (np.max(diff) < 1e-5): raise AssertionError() + if not (np.max(diff) < 1e-5): + raise AssertionError() model.add('c1_amp', 1, vary=True) model.add('c1_xo', 5, vary=True) @@ -109,7 +120,8 @@ def test_emp_vs_ana_jacobian(): emp_Jij = fitting.emp_jacobian(model, x, y) ana_Jij = fitting.jacobian(model, x, y) diff = np.abs(ana_Jij - emp_Jij) - if not (np.max(diff) < 1e-3): raise AssertionError() + if not (np.max(diff) < 1e-3): + raise AssertionError() def test_emp_vs_ana_hessian(): @@ -120,7 +132,8 @@ def test_emp_vs_ana_hessian(): emp_Hij = fitting.emp_hessian(model, x, y) ana_Hij = fitting.hessian(model, x, y) diff = np.abs(ana_Hij - emp_Hij) - if not (np.max(diff) < 1e-5): raise AssertionError() + if not (np.max(diff) < 1e-5): + raise AssertionError() model.add('c1_amp', 1, vary=True) model.add('c1_xo', 5, vary=True) @@ -133,17 +146,20 @@ def test_emp_vs_ana_hessian(): emp_Hij = fitting.emp_hessian(model, x, y) ana_Hij = fitting.hessian(model, x, y) diff = np.abs(ana_Hij - emp_Hij) - if not (np.max(diff) < 1): raise AssertionError() + if not (np.max(diff) < 1): + raise AssertionError() def test_make_ita(): """Test make_ita""" noise = np.random.random((10, 10)) ita = fitting.make_ita(noise) - if not (ita.shape == (100, 100)): raise AssertionError() + if not (ita.shape == (100, 100)): + raise AssertionError() noise *= np.nan ita = fitting.make_ita(noise) - if not (len(ita) == 0): raise AssertionError() + if not (len(ita) == 0): + raise AssertionError() def test_RB_bias(): @@ -151,7 +167,8 @@ def test_RB_bias(): data = np.random.random((4, 4)) model = make_model() bias = fitting.RB_bias(data, model) - if not (len(bias) == 3): raise AssertionError() + if not (len(bias) == 3): + raise AssertionError() def test_bias_correct(): @@ -172,15 +189,24 @@ def test_condon_errs(): source.local_rms = 0.1 source.peak_flux = 1 source.int_flux = 1 + fitting.condon_errors(source, None) - if not (source.err_a == 0): raise AssertionError() + if not (source.err_a == 0): + raise AssertionError( + "err_a should be zero but is {0}".format(source.err_a)) + fitting.condon_errors(source, theta_n=8.) - if not (source.err_a > 0): raise AssertionError() + if not (source.err_a > 0): + raise AssertionError("err_a should be non-zero") + # test that we can get a PA error source.a = 20 fitting.condon_errors(source, None) + if source.err_pa < 0: + raise AssertionError("err_pa<0") if source.err_pa < 0: raise AssertionError() + if __name__ == "__main__": # introspect and run all the functions starting with 'test' for f in dir(): From fc052ad33057b6a13a3a68c826c6cfc181d3d67e Mon Sep 17 00:00:00 2001 From: PaulHancock Date: Mon, 15 Nov 2021 14:06:46 +0800 Subject: [PATCH 2/3] detect and fix bug Bug when using the --condon option #156 --- AegeanTools/fitting.py | 5 ++--- tests/test_fitting.py | 9 +++++++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/AegeanTools/fitting.py b/AegeanTools/fitting.py index ac7e29e9..41523563 100644 --- a/AegeanTools/fitting.py +++ b/AegeanTools/fitting.py @@ -883,9 +883,8 @@ def condon_errors(source, theta_n, psf=None): minor = source.b / 3600. # degrees phi = np.radians(source.pa) # radians if psf is not None: - beam = psf.get_beam(source.ra, source.dec) - if beam is not None: - theta_n = np.sqrt(beam.a * beam.b) + a, b, _ = psf.get_psf_sky2sky(source.ra, source.dec) + theta_n = np.sqrt(a * b) if theta_n is None: source.err_a = source.err_b = source.err_peak_flux = source.err_pa = source.err_int_flux = 0.0 diff --git a/tests/test_fitting.py b/tests/test_fitting.py index 5f71aca5..875033f8 100755 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -7,7 +7,7 @@ __author__ = 'Paul Hancock' -from AegeanTools import fitting, models +from AegeanTools import fitting, models, wcs_helpers import lmfit import numpy as np @@ -204,7 +204,12 @@ def test_condon_errs(): fitting.condon_errors(source, None) if source.err_pa < 0: raise AssertionError("err_pa<0") - if source.err_pa < 0: raise AssertionError() + + psf = wcs_helpers.WCSHelper.from_file('tests/test_files/1904-66_SIN.fits') + try: + fitting.condon_errors(source, theta_n=8., psf=psf) + except AttributeError as e: + raise AssertionError("condon_errors failed with psf") if __name__ == "__main__": From a02c2cd2c65dacbc49aea24f6e52a29d7b1a7720 Mon Sep 17 00:00:00 2001 From: PaulHancock Date: Mon, 15 Nov 2021 14:09:04 +0800 Subject: [PATCH 3/3] update date --- AegeanTools/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/AegeanTools/__init__.py b/AegeanTools/__init__.py index 365a02b2..63d5155a 100755 --- a/AegeanTools/__init__.py +++ b/AegeanTools/__init__.py @@ -11,7 +11,7 @@ """ __author__ = 'Paul Hancock' __version__ = '2.2.5' -__date__ = '2021-11-10' +__date__ = '2021-11-15' __citation__ = """ % If your work makes use of AegeanTools please cite the following papers as appropriate: