<commit_before>"""
Some generic utility routines for number handling and
calculating (specific) variances
"""

import itertools
import logging

import numpy

from sourcefinder import extract
from sourcefinder import stats
from sourcefinder import utils
from sourcefinder.utility import containers
from sourcefinder.utility.memoize import Memoize

import psutil
import time
import dask.array as da

try:
    import ndimage
except ImportError:
    from scipy import ndimage
    
def timeit(method):
    def timed(*args, **kw):
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()
        if 'log_time' in kw:
            name = kw.get('log_name', method._name_.upper())
            kw['log_time'][name] = int((te - ts) * 1000)
        else:
            print('{0}  {1:2.2f} ms'.format(method.__name__, (te - ts) * 1000))
        return result
    return timed

logger = logging.getLogger(__name__)

#
# Hard-coded configuration parameters; not user settable.
#
INTERPOLATE_ORDER = 1  # Spline order for grid interpolation
MEDIAN_FILTER = 0  # If non-zero, apply a median filter of size
# MEDIAN_FILTER to the background and RMS grids prior
# to interpolating.
MF_THRESHOLD = 0  # If MEDIAN_FILTER is non-zero, only use the filtered
# grid when the (absolute) difference between the raw
# and filtered grids is larger than MF_THRESHOLD.
DEBLEND_MINCONT = 0.005  # Min. fraction of island flux in deblended subisland
STRUCTURING_ELEMENT = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]  # Island connectiivty

class ImageData(object):
    """Encapsulates an image in terms of a numpy array + meta/headerdata.

    This is your primary contact point for interaction with images: it icludes
    facilities for source extraction and measurement, etc.
    """

    def __init__(self, data, beam, wcs, margin=0, radius=0, back_size_x=32,
                 back_size_y=32, residuals=True
                 ):
        """Sets up an ImageData object.

        *Args:*
          - data (2D numpy.ndarray): actual image data
          - wcs (utility.coordinates.wcs): world coordinate system
            specification
          - beam (3-tuple): beam shape specification as
            (semimajor, semiminor, theta)

        """

        # Do data, wcs and beam need deepcopy?
        # Probably not (memory overhead, in particular for data),
        # but then the user shouldn't change them outside ImageData in the
        # mean time
        self.rawdata = data  # a 2D numpy array
        self.wcs = wcs  # a utility.coordinates.wcs instance
        self.beam = beam  # tuple of (semimaj, semimin, theta)
        self.clip = {}
        self.labels = {}
        self.freq_low = 1
        self.freq_high = 1

        self.back_size_x = back_size_x
        self.back_size_y = back_size_y
        self.margin = margin
        self.radius = radius
        self.residuals = residuals

    ###########################################################################
    #                                                                         #
    # Properties and attributes.                                              #
    #                                                                         #
    # Properties are attributes managed by methods; rather than calling the   #
    # method directly, the attribute automatically invokes it. We can use     #
    # this to do cunning transparent caching ("memoizing") etc; see the       #
    # Memoize class.                                                          #
    #                                                                         #
    # clearcache() clears all the memoized data, which can get quite large.   #
    # It may be wise to call this, for example, in an exception handler       #
    # dealing with MemoryErrors.                                              #
    #                                                                         #
    ###########################################################################
    @Memoize
    def _grids(self):
        """Gridded RMS and background data for interpolating"""
        return self.__grids()

    grids = property(fget=_grids, fdel=_grids.delete)

    @Memoize
    def _backmap(self):
        """Background map"""
        if not hasattr(self, "_user_backmap"):
            return self._interpolate(self.grids['bg'])
        else:
            return self._user_backmap

    def _set_backmap(self, bgmap):
        self._user_backmap = bgmap
        del (self.backmap)
        del (self.data_bgsubbed)

    backmap = property(fget=_backmap, fdel=_backmap.delete, fset=_set_backmap)

    @Memoize
    def _get_rm(self):
        """RMS map"""
        if not hasattr(self, "_user_noisemap"):
            return self._interpolate(self.grids['rms'], roundup=True)
        else:
            return self._user_noisemap

    def _set_rm(self, noisemap):
        self._user_noisemap = noisemap
        del (self.rmsmap)

    rmsmap = property(fget=_get_rm, fdel=_get_rm.delete, fset=_set_rm)

    @Memoize
    def _get_data(self):
        """Masked image data"""
        # We will ignore all the data which is masked for the rest of the
        # sourcefinding process. We build up the mask by stacking ("or-ing
        # together") a number of different effects:
        #
        # * A margin from the edge of the image;
        # * Any data outside a given radius from the centre of the image;
        # * Data which is "obviously" bad (equal to 0 or NaN).
        mask = numpy.zeros((self.xdim, self.ydim))
        if self.margin:
            margin_mask = numpy.ones((self.xdim, self.ydim))
            margin_mask[self.margin:-self.margin, self.margin:-self.margin] = 0
            mask = numpy.logical_or(mask, margin_mask)
        if self.radius:
            radius_mask = utils.circular_mask(self.xdim, self.ydim, self.radius)
            mask = numpy.logical_or(mask, radius_mask)
        mask = numpy.logical_or(mask, numpy.isnan(self.rawdata))
        return numpy.ma.array(self.rawdata, mask=mask)

    data = property(fget=_get_data, fdel=_get_data.delete)

    @Memoize
    def _get_data_bgsubbed(self):
        """Background subtracted masked image data"""
        return self.data - self.backmap

    data_bgsubbed = property(fget=_get_data_bgsubbed,
                             fdel=_get_data_bgsubbed.delete)

    @property
    def xdim(self):
        """X pixel dimension of (unmasked) data"""
        return self.rawdata.shape[0]

    @property
    def ydim(self):
        """Y pixel dimension of (unmasked) data"""
        return self.rawdata.shape[1]

    @property
    def pixmax(self):
        """Maximum pixel value (pre-background subtraction)"""
        return self.data.max()

    @property
    def pixmin(self):
        """Minimum pixel value (pre-background subtraction)"""
        return self.data.min()

    def clearcache(self):
        """Zap any calculated data stored in this object.

        Clear the background and rms maps, labels, clip, and any locally held
        data. All of these can be reconstructed from the data accessor.

        Note that this *must* be run to pick up any new settings.
        """
        self.labels.clear()
        self.clip.clear()
        del (self.backmap)
        del (self.rmsmap)
        del (self.data)
        del (self.data_bgsubbed)
        del (self.grids)
        if hasattr(self, 'residuals_from_gauss_fitting'):
            del (self.residuals_from_gauss_fitting)
        if hasattr(self, 'residuals_from_deblending'):
            del (self.residuals_from_deblending)

    ###########################################################################
    #                                                                         #
    # General purpose image handling.                                         #
    #                                                                         #
    # Routines for saving and trimming data, and calculating background/RMS   #
    # maps (in conjuntion with the properties above).                         #
    #                                                                         #
    ###########################################################################

    # Private "support" methods
    @timeit
    def __grids(self):
        """Calculate background and RMS grids of this image.

        These grids can be interpolated up to make maps of the original image
        dimensions: see _interpolate().

        This is called automatically when ImageData.backmap,
        ImageData.rmsmap or ImageData.fdrmap is first accessed.
        """

        # there's no point in working with the whole of the data array
        # if it's masked.
        useful_chunk = ndimage.find_objects(numpy.where(self.data.mask, 0, 1))
        assert (len(useful_chunk) == 1)
        useful_data = da.from_array(self.data[useful_chunk[0]].data, chunks=(self.back_size_x, self.back_size_y))

        mode_and_rms = useful_data.map_blocks(self.compute_mode_and_rms_of_subimages, dtype=numpy.complex64,
                                              chunks=(1, 1)).compute()

        # See also similar comment below. This solution was chosen because map_blocks does not seem to be able to
        # output multiple arrays. One can however output to a complex array and take real and imaginary
        # parts afterwards. Not a very clean solution, I admit.
        mode_grid = mode_and_rms.real
        rms_grid = mode_and_rms.imag

        rms_grid = numpy.ma.array(
            rms_grid, mask=numpy.where(rms_grid == 0, 1, 0))
        # A rms of zero is not physical, since any instrument has system noise, so I use that as criterion<commit_msg>Seems more apprpriate to use the class name (ImageData) instead of self, the class instantiation, but perhaps the former applies to class methods instead of static methods.<commit_after>"""
Some generic utility routines for number handling and
calculating (specific) variances
"""

import itertools
import logging

import numpy

from sourcefinder import extract
from sourcefinder import stats
from sourcefinder import utils
from sourcefinder.utility import containers
from sourcefinder.utility.memoize import Memoize

import psutil
import time
import dask.array as da

try:
    import ndimage
except ImportError:
    from scipy import ndimage
    
def timeit(method):
    def timed(*args, **kw):
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()
        if 'log_time' in kw:
            name = kw.get('log_name', method._name_.upper())
            kw['log_time'][name] = int((te - ts) * 1000)
        else:
            print('{0}  {1:2.2f} ms'.format(method.__name__, (te - ts) * 1000))
        return result
    return timed

logger = logging.getLogger(__name__)

#
# Hard-coded configuration parameters; not user settable.
#
INTERPOLATE_ORDER = 1  # Spline order for grid interpolation
MEDIAN_FILTER = 0  # If non-zero, apply a median filter of size
# MEDIAN_FILTER to the background and RMS grids prior
# to interpolating.
MF_THRESHOLD = 0  # If MEDIAN_FILTER is non-zero, only use the filtered
# grid when the (absolute) difference between the raw
# and filtered grids is larger than MF_THRESHOLD.
DEBLEND_MINCONT = 0.005  # Min. fraction of island flux in deblended subisland
STRUCTURING_ELEMENT = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]  # Island connectiivty

class ImageData(object):
    """Encapsulates an image in terms of a numpy array + meta/headerdata.

    This is your primary contact point for interaction with images: it icludes
    facilities for source extraction and measurement, etc.
    """

    def __init__(self, data, beam, wcs, margin=0, radius=0, back_size_x=32,
                 back_size_y=32, residuals=True
                 ):
        """Sets up an ImageData object.

        *Args:*
          - data (2D numpy.ndarray): actual image data
          - wcs (utility.coordinates.wcs): world coordinate system
            specification
          - beam (3-tuple): beam shape specification as
            (semimajor, semiminor, theta)

        """

        # Do data, wcs and beam need deepcopy?
        # Probably not (memory overhead, in particular for data),
        # but then the user shouldn't change them outside ImageData in the
        # mean time
        self.rawdata = data  # a 2D numpy array
        self.wcs = wcs  # a utility.coordinates.wcs instance
        self.beam = beam  # tuple of (semimaj, semimin, theta)
        self.clip = {}
        self.labels = {}
        self.freq_low = 1
        self.freq_high = 1

        self.back_size_x = back_size_x
        self.back_size_y = back_size_y
        self.margin = margin
        self.radius = radius
        self.residuals = residuals

    ###########################################################################
    #                                                                         #
    # Properties and attributes.                                              #
    #                                                                         #
    # Properties are attributes managed by methods; rather than calling the   #
    # method directly, the attribute automatically invokes it. We can use     #
    # this to do cunning transparent caching ("memoizing") etc; see the       #
    # Memoize class.                                                          #
    #                                                                         #
    # clearcache() clears all the memoized data, which can get quite large.   #
    # It may be wise to call this, for example, in an exception handler       #
    # dealing with MemoryErrors.                                              #
    #                                                                         #
    ###########################################################################
    @Memoize
    def _grids(self):
        """Gridded RMS and background data for interpolating"""
        return self.__grids()

    grids = property(fget=_grids, fdel=_grids.delete)

    @Memoize
    def _backmap(self):
        """Background map"""
        if not hasattr(self, "_user_backmap"):
            return self._interpolate(self.grids['bg'])
        else:
            return self._user_backmap

    def _set_backmap(self, bgmap):
        self._user_backmap = bgmap
        del (self.backmap)
        del (self.data_bgsubbed)

    backmap = property(fget=_backmap, fdel=_backmap.delete, fset=_set_backmap)

    @Memoize
    def _get_rm(self):
        """RMS map"""
        if not hasattr(self, "_user_noisemap"):
            return self._interpolate(self.grids['rms'], roundup=True)
        else:
            return self._user_noisemap

    def _set_rm(self, noisemap):
        self._user_noisemap = noisemap
        del (self.rmsmap)

    rmsmap = property(fget=_get_rm, fdel=_get_rm.delete, fset=_set_rm)

    @Memoize
    def _get_data(self):
        """Masked image data"""
        # We will ignore all the data which is masked for the rest of the
        # sourcefinding process. We build up the mask by stacking ("or-ing
        # together") a number of different effects:
        #
        # * A margin from the edge of the image;
        # * Any data outside a given radius from the centre of the image;
        # * Data which is "obviously" bad (equal to 0 or NaN).
        mask = numpy.zeros((self.xdim, self.ydim))
        if self.margin:
            margin_mask = numpy.ones((self.xdim, self.ydim))
            margin_mask[self.margin:-self.margin, self.margin:-self.margin] = 0
            mask = numpy.logical_or(mask, margin_mask)
        if self.radius:
            radius_mask = utils.circular_mask(self.xdim, self.ydim, self.radius)
            mask = numpy.logical_or(mask, radius_mask)
        mask = numpy.logical_or(mask, numpy.isnan(self.rawdata))
        return numpy.ma.array(self.rawdata, mask=mask)

    data = property(fget=_get_data, fdel=_get_data.delete)

    @Memoize
    def _get_data_bgsubbed(self):
        """Background subtracted masked image data"""
        return self.data - self.backmap

    data_bgsubbed = property(fget=_get_data_bgsubbed,
                             fdel=_get_data_bgsubbed.delete)

    @property
    def xdim(self):
        """X pixel dimension of (unmasked) data"""
        return self.rawdata.shape[0]

    @property
    def ydim(self):
        """Y pixel dimension of (unmasked) data"""
        return self.rawdata.shape[1]

    @property
    def pixmax(self):
        """Maximum pixel value (pre-background subtraction)"""
        return self.data.max()

    @property
    def pixmin(self):
        """Minimum pixel value (pre-background subtraction)"""
        return self.data.min()

    def clearcache(self):
        """Zap any calculated data stored in this object.

        Clear the background and rms maps, labels, clip, and any locally held
        data. All of these can be reconstructed from the data accessor.

        Note that this *must* be run to pick up any new settings.
        """
        self.labels.clear()
        self.clip.clear()
        del (self.backmap)
        del (self.rmsmap)
        del (self.data)
        del (self.data_bgsubbed)
        del (self.grids)
        if hasattr(self, 'residuals_from_gauss_fitting'):
            del (self.residuals_from_gauss_fitting)
        if hasattr(self, 'residuals_from_deblending'):
            del (self.residuals_from_deblending)

    ###########################################################################
    #                                                                         #
    # General purpose image handling.                                         #
    #                                                                         #
    # Routines for saving and trimming data, and calculating background/RMS   #
    # maps (in conjuntion with the properties above).                         #
    #                                                                         #
    ###########################################################################

    # Private "support" methods
    @timeit
    def __grids(self):
        """Calculate background and RMS grids of this image.

        These grids can be interpolated up to make maps of the original image
        dimensions: see _interpolate().

        This is called automatically when ImageData.backmap,
        ImageData.rmsmap or ImageData.fdrmap is first accessed.
        """

        # there's no point in working with the whole of the data array
        # if it's masked.
        useful_chunk = ndimage.find_objects(numpy.where(self.data.mask, 0, 1))
        assert (len(useful_chunk) == 1)
        useful_data = da.from_array(self.data[useful_chunk[0]].data, chunks=(self.back_size_x, self.back_size_y))

        mode_and_rms = useful_data.map_blocks(ImageData.compute_mode_and_rms_of_subimages, dtype=numpy.complex64,
                                              chunks=(1, 1)).compute()

        # See also similar comment below. This solution was chosen because map_blocks does not seem to be able to
        # output multiple arrays. One can however output to a complex array and take real and imaginary
        # parts afterwards. Not a very clean solution, I admit.
        mode_grid = mode_and_rms.real
        rms_grid = mode_and_rms.imag

        rms_grid = numpy.ma.array(
            rms_grid, mask=numpy.where(rms_grid == 0, 1, 0))
        # A rms of zero is not physical, since any instrument has system noise, so I use that as criterion
{'fraction_non_alphanumeric': 0.11085556577736891, 'fraction_numerical': 0.005519779208831647, 'mean_word_length': 2.0056835637480797, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 60, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3703869', 'n_tokens_mistral': 5438, 'n_tokens_neox': 5215, 'n_words': 1661}