<commit_before>"""
Documentation.
"""

import numpy
import scipy.signal
import scipy.integrate
from abc import abstractmethod

import ugali.utils.plotting
import ugali.utils.projector

############################################################

class Kernel:
    """
    Normalized spatial kernel to describe candidate clusters.
    """
    
    def __init__(self, config):

        self.config = config

        if self.config.params['kernel']['type'].lower() == 'disk':
            self.kernel, self.bins_x, self.bins_y = disk(self.config.params['kernel']['params'],
                                                         self.config.params['coords']['pixel_size'])
        else:
            print 'WARNING: kernel type %s not recognized'%(self.config.params['kernel']['type'])

    def plot(self):
        ugali.utils.plotting.twoDimensionalHistogram('kernel', 'x (deg)', 'y (deg)',
                                                     self.kernel, self.bins_x, self.bins_y,
                                                     lim_x = [self.bins_x[0], self.bins_x[-1]],
                                                     lim_y = [self.bins_y[0], self.bins_y[-1]])
    
    #def convolve(self, map):
    #    map_convolve = scipy.signal.convolve(map, self.kernel, mode='same')
    #    return map_convolve()
        

############################################################

def disk(params, pixel_size):

    disk_radius = params[0]

    disk_radius_pixels = numpy.ceil(disk_radius / pixel_size) # Disk radius in projected map pixels

    bins_x = numpy.linspace(-pixel_size * (disk_radius_pixels + 0.5),
                            pixel_size * (disk_radius_pixels + 0.5),
                            2 * (int(disk_radius_pixels) + 1))
    bins_y = bins_x
    centers_x = bins_x[0: -1] + (0.5 * pixel_size)
    centers_y = centers_x
    
    mesh_x, mesh_y = numpy.meshgrid(centers_x, centers_y)

    r = numpy.sqrt(mesh_x**2 + mesh_y**2)

    kernel = r < disk_radius

    kernel = kernel.astype(float)

    kernel /= float(numpy.sum(kernel))

    return kernel, bins_x, bins_y

############################################################

class RadialKernel(object):
    """
    Base class for radially symmetric kernel profiles.
    Each subclass must implement:
        pdf       : value of the pdf as a function of angular radius
        extension : parameter for controlling extension
        edge      : parameter for defining effective edge of the kernel
    """
    def __init__(self, lon, lat, *args):
        self.name = self.__class__.__name__

        self.setCenter(lon,lat)
        self.setParams(*args)

    def setParams(self, *args):
        self.params = args
        self.norm = 1.0
        #self.norm /= self.integrate()

    def __call__(self, radius):
        return self.pdf(radius)

    @abstractmethod
    def extension(self):
        pass

    @abstractmethod
    def edge(self):
        pass

    @abstractmethod
    def pdf(self, radius):
        """
        Evaluate the PDF (deg^2) at a given radius (deg).
        """
        pass

    def integrate(self, r_min=0, r_max=numpy.inf):
        """
        Calculate the 2D integral of the PDF.
        """
        if r_min < 0: raise Exception('r_min must be >= 0')
        r_max = r_max if r_max < self.edge() else self.edge()
        r_min,r_max = numpy.radians([r_min,r_max])
        integrand = lambda r: self.pdf(numpy.degrees(r)) * 2*numpy.pi*r
        return scipy.integrate.quad(integrand, r_min, r_max, full_output=True, epsabs=0)[0]

    # For back-compatibility
    def integratePDF(self, r_min, r_max, steps=1e4, prenormalized=True): 
        return self.integrate(r_min, r_max)

    # For back-compatibility
    def surfaceIntensity(self, radius): 
        return self.pdf(radius)

    def setCenter(self, lon, lat):
        self.lon = lon
        self.lat = lat
        self.projector = ugali.utils.projector.Projector(self.lon, self.lat)

    def sample_radius(self, n):
        """
        Sample the radial distribution (deg) from the 2D stellar density.
        """
        radius = numpy.linspace(0, self.edge(), 1.e5)
        rpdf = self.pdf(r) * numpy.sin(numpy.radians(radius))
        rcdf = numpy.cumsum(rpdf)
        rcdf /= rcdf[-1]
        fn = scipy.interpolate.interp1d(rcdf, range(0, len(rcdf)), bounds_error=False, fill_value=-1)
        index = numpy.floor(fn(numpy.random.rand(n))).astype(int)
        return radius[index]

    def sample_lonlat(self, n):
        """
        Sample 2D distribution of points in lon, lat
        Careful, doesn't respect sky projections
        """
        radius = self.sample_radius(n)<commit_msg>Fixed integration to work in degrees<commit_after>"""
Documentation.
"""

import numpy
import scipy.signal
import scipy.integrate
from abc import abstractmethod

import ugali.utils.plotting
import ugali.utils.projector

############################################################

class Kernel:
    """
    Normalized spatial kernel to describe candidate clusters.
    """
    
    def __init__(self, config):

        self.config = config

        if self.config.params['kernel']['type'].lower() == 'disk':
            self.kernel, self.bins_x, self.bins_y = disk(self.config.params['kernel']['params'],
                                                         self.config.params['coords']['pixel_size'])
        else:
            print 'WARNING: kernel type %s not recognized'%(self.config.params['kernel']['type'])

    def plot(self):
        ugali.utils.plotting.twoDimensionalHistogram('kernel', 'x (deg)', 'y (deg)',
                                                     self.kernel, self.bins_x, self.bins_y,
                                                     lim_x = [self.bins_x[0], self.bins_x[-1]],
                                                     lim_y = [self.bins_y[0], self.bins_y[-1]])
    
    #def convolve(self, map):
    #    map_convolve = scipy.signal.convolve(map, self.kernel, mode='same')
    #    return map_convolve()
        

############################################################

def disk(params, pixel_size):

    disk_radius = params[0]

    disk_radius_pixels = numpy.ceil(disk_radius / pixel_size) # Disk radius in projected map pixels

    bins_x = numpy.linspace(-pixel_size * (disk_radius_pixels + 0.5),
                            pixel_size * (disk_radius_pixels + 0.5),
                            2 * (int(disk_radius_pixels) + 1))
    bins_y = bins_x
    centers_x = bins_x[0: -1] + (0.5 * pixel_size)
    centers_y = centers_x
    
    mesh_x, mesh_y = numpy.meshgrid(centers_x, centers_y)

    r = numpy.sqrt(mesh_x**2 + mesh_y**2)

    kernel = r < disk_radius

    kernel = kernel.astype(float)

    kernel /= float(numpy.sum(kernel))

    return kernel, bins_x, bins_y

############################################################

class RadialKernel(object):
    """
    Base class for radially symmetric kernel profiles.
    Each subclass must implement:
        pdf       : value of the pdf as a function of angular radius
        extension : parameter for controlling extension
        edge      : parameter for defining effective edge of the kernel
    """
    def __init__(self, lon, lat, *args):
        self.name = self.__class__.__name__

        self.setCenter(lon,lat)
        self.setParams(*args)

    def setParams(self, *args):
        self.params = args
        self.norm = 1.0
        self.norm /= self.integrate()

    def __call__(self, radius):
        return self.pdf(radius)

    @abstractmethod
    def extension(self):
        pass

    @abstractmethod
    def edge(self):
        pass

    @abstractmethod
    def pdf(self, radius):
        """
        Evaluate the surface brightness profile, I(R), at a given radius (deg).
        The surface brightness profiles should be defined in the manner of
        second convention of Binney and Tremain, 2008 (Box 2.1)
        radius : radius to evaluate the local surface brightness
        return : surface brightness (1/deg^2)
        """
        pass

    # For back-compatibility
    def surfaceIntensity(self, radius): 
        return self.pdf(radius)

    def integrate(self, r_min=0, r_max=numpy.inf):
        """
        Calculate the 2D integral of the surface brightness profile (i.e, the flux) 
        between r_min and r_max.
        r_min : minimum integration radius
        r_max : maximum integration radius
        return : Solid angle integral 
        """
        if r_min < 0: raise Exception('r_min must be >= 0')
        r_max = r_max if r_max < self.edge() else self.edge()
        r_min,r_max = numpy.radians([r_min,r_max])
        integrand = lambda r: self.pdf(r) * 2*numpy.pi * r
        return scipy.integrate.quad(integrand, r_min, r_max, full_output=True, epsabs=0)[0]

    # For back-compatibility
    def integratePDF(self, r_min, r_max, steps=1e4, prenormalized=True): 
        return self.integrate(r_min, r_max)

    def setCenter(self, lon, lat):
        self.lon = lon
        self.lat = lat
        self.projector = ugali.utils.projector.Projector(self.lon, self.lat)

    def sample_radius(self, n):
        """
        Sample the radial distribution (deg) from the 2D stellar density.
        """
        radius = numpy.linspace(0, self.edge(), 1.e5)
        rpdf = self.pdf(r) * numpy.sin(numpy.radians(radius))
        rcdf = numpy.cumsum(rpdf)
        rcdf /= rcdf[-1]
        fn = scipy.interpolate.interp1d(rcdf, range(0, len(rcdf)), bounds_error=False, fill_value=-1)
        index = numpy.floor(fn(numpy.random.rand(n))).astype(int)
        return radius[index]

    def sample_lonlat(self, n):
        """
        Sample 2D distribution of points in lon, lat
        Careful, doesn't respect sky projections
        """
        radius = self.sample_radius(n)
{'fraction_non_alphanumeric': 0.1446938775510204, 'fraction_numerical': 0.00846938775510204, 'mean_word_length': 2.2955615332885, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3663937', 'n_tokens_mistral': 3023, 'n_tokens_neox': 2882, 'n_words': 710}