remove unused code related to bias correction

AegeanTools/fitting.py

@@ -691,188 +691,6 @@ def emp_hessian(pars, x, y):
     return matrix
 
 
-def nan_acf(noise):
-    """
-    Calculate the autocorrelation function of the noise
-    where the noise is a 2d array that may contain nans
-
-
-    Parameters
-    ----------
-    noise : 2d-array
-        Noise image.
-
-    Returns
-    -------
-    acf : 2d-array
-        The ACF.
-    """
-    corr = np.zeros(noise.shape)
-    ix, jx = noise.shape
-    for i in range(ix):
-        si_min = slice(i, None, None)
-        si_max = slice(None, ix - i, None)
-        for j in range(jx):
-            sj_min = slice(j, None, None)
-            sj_max = slice(None, jx - j, None)
-            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])
-    # return the normalised acf
-    return corr / np.nanmax(corr)
-
-
-def make_ita(noise, acf=None):
-    """
-    Create the matrix ita of the noise where the noise may be a masked array
-    where ita(x,y) is the correlation between pixel pairs that have the same
-    separation as x and y.
-
-    Parameters
-    ----------
-    noise : 2d-array
-        The noise image
-
-    acf : 2d-array
-        The autocorrelation matrix. (None = calculate from data).
-        Default = None.
-
-    Returns
-    -------
-    ita : 2d-array
-        The matrix ita
-    """
-    if acf is None:
-        acf = nan_acf(noise)
-    # s should be the number of non-masked pixels
-    s = np.count_nonzero(np.isfinite(noise))
-    # the indices of the non-masked pixels
-    xm, ym = np.where(np.isfinite(noise))
-    ita = np.zeros((s, s))
-    # iterate over the pixels
-    for i, (x1, y1) in enumerate(zip(xm, ym)):
-        for j, (x2, y2) in enumerate(zip(xm, ym)):
-            k = abs(x1 - x2)
-            l = abs(y1 - y2)
-            ita[i, j] = acf[k, l]
-    return ita
-
-
-def RB_bias(data, pars, ita=None, acf=None):
-    """
-    Calculate the expected bias on each of the parameters in the model pars.
-    Only parameters that are allowed to vary will have a bias.
-    Calculation follows the description of Refrieger & Brown 1998 (cite).
-
-
-    Parameters
-    ----------
-    data : 2d-array
-        data that was fit
-
-    pars : lmfit.Parameters
-        The model
-
-    ita : 2d-array
-        The ita matrix (optional).
-
-    acf : 2d-array
-        The acf for the data.
-
-    Returns
-    -------
-    bias : array
-        The bias on each of the parameters
-    """
-    logger.info("data {0}".format(data.shape))
-    nparams = np.sum([pars[k].vary for k in pars.keys() if k != "components"])
-    # masked pixels
-    xm, ym = np.where(np.isfinite(data))
-    # 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
-    )
-
-    h = hessian(pars, x, y)
-    # mask the hessian to be AxAxN array
-    h = h[:, :, xm, ym]
-    Hij = np.einsum("ik,jk", j, j)
-    Dij = np.linalg.inv(Hij)
-    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_2 = -1.0 / 2 * np.einsum("rkj,ir,km,jn", Bijk, Dij, Dij, Dij)
-    Cimn = Cimn_1 + Cimn_2
-
-    if ita is None:
-        # N is the noise (data-model)
-        N = data - ntwodgaussian_lmfit(pars)(x, y)
-        if acf is None:
-            acf = nan_acf(N)
-        ita = make_ita(N, acf=acf)
-        logger.info("acf.shape {0}".format(acf.shape))
-        logger.info("acf[0] {0}".format(acf[0]))
-        logger.info("ita.shape {0}".format(ita.shape))
-        logger.info("ita[0] {0}".format(ita[0]))
-
-    # Included for completeness but not required
-
-    # now mask/ravel the noise
-    # N = N[np.isfinite(N)].ravel()
-    # Pi = np.einsum('ip,p', j, N)
-    # Qij = np.einsum('ijp,p', h, N)
-
-    Vij = np.einsum("ip,jq,pq", j, j, ita)
-    Uijk = np.einsum("ip,jkq,pq", j, h, ita)
-
-    bias_1 = np.einsum("imn, mn", Cimn, Vij)
-    bias_2 = np.einsum("ilkm, mlk", Eilkm, Uijk)
-    bias = bias_1 + bias_2
-    logger.info("bias {0}".format(bias))
-    return bias
-
-
-def bias_correct(params, data, acf=None):
-    """
-    Calculate and apply a bias correction to the given fit parameters
-
-
-    Parameters
-    ----------
-    params : lmfit.Parameters
-        The model parameters. These will be modified.
-
-    data : 2d-array
-        The data which was used in the fitting
-
-    acf : 2d-array
-        ACF of the data. Default = None.
-
-    Returns
-    -------
-    None
-
-    See Also
-    --------
-    :func:`AegeanTools.fitting.RB_bias`
-    """
-    bias = RB_bias(data, params, acf=acf)
-    i = 0
-    for p in params:
-        if "theta" in p:
-            continue
-        if params[p].vary:
-            params[p].value -= bias[i]
-            i += 1
-    return
-
-
 def errors(source, model, wcshelper):
     """
     Convert pixel based errors into sky coord errors

AegeanTools/source_finder.py

@@ -24,7 +24,6 @@
 from .fitting import (
     Bmatrix,
     Cmatrix,
-    bias_correct,
     covar_errors,
     do_lmfit,
     do_lmfit_3D,
@@ -192,9 +191,6 @@ class SourceFinder(object):
     docov : bool
         If True, then fitting will be done using the covariance matrix.
 
-    dobais: bool
-        Not yet implemented. Default = False
-
     sources : [ComponentSource|IslandSource|SimpleSource]
       List of sources that have been found/measured.
     """
@@ -212,7 +208,6 @@ def __init__(self, **kwargs):
         self.wcshelper = None
         self.blank = False
         self.docov = True
-        self.dobias = False
         self.cube_index = 0
 
         self.sources = []
@@ -1049,9 +1044,6 @@ def load_globals(
         self.blank = blank
         self.docov = docov
 
-        # Default to false until I can verify that this is working
-        self.dobias = False
-
         # check if the WCS is galactic
         if "lon" in self.header["CTYPE1"].lower():
             logger.info("Galactic coordinates detected and noted")
@@ -2121,20 +2113,6 @@ def _fit_island(self, island_data):
                 # get the real (sky) parameter errors
                 model = covar_errors(result.params, idata, errs=errs, B=B, C=C)
 
-                if self.dobias and self.docov:
-                    x, y = np.indices(idata.shape)
-                    acf = elliptical_gaussian(
-                        x,
-                        y,
-                        1,
-                        0,
-                        0,
-                        pixbeam.a * FWHM2CC * fac,
-                        pixbeam.b * FWHM2CC * fac,
-                        pixbeam.pa,
-                    )
-                    bias_correct(model, idata, acf=acf * errs**2)
-
                 if not result.success:
                     is_flag |= flags.FITERR