Gabriel-p/KDE_mult_overdens.py

## KDE_mult_overdens.py

import numpy as np
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.filters import gaussian_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.offsetbox as offsetbox
# from mpl_toolkits.axes_grid1 import make_axes_locatable


def get_coords():
    '''
    Return coordinates from file.
    '''
    # Each sub-list in file_data is a row of the file.
    # file_data = np.loadtxt("TR24_clean.dat")
    file_data = np.loadtxt("CLUSTER.DAT")

    # Extract coordinates and zip them into two lists.
    ra, dec = zip(*file_data)[1], zip(*file_data)[2]

    # ra_rang, dec_rang = max(ra) - min(ra), max(dec) - min(dec)
    # print 'RA range: ', ra_rang
    # print 'DEC range: ', dec_rang

    return ra, dec


def get_2d_histo(x_data, y_data, N_bins, std_dev):
    '''
    Return 2D histogram for the positional data.
    '''
    xmin, xmax = min(x_data), max(x_data)
    ymin, ymax = min(y_data), max(y_data)

    # Calculate the number of bins used.
    x_rang, y_rang = (xmax - xmin), (ymax - ymin)

    # Bin width to create the 2D histogram.
    bin_width = min(x_rang, y_rang) / N_bins

    # Number of bins in x,y given the bin width.
    binsxy = [int(x_rang / bin_width), int(y_rang / bin_width)]

    # hist_2d is the 2D histogram, *edges store the edges of the bins.
    hist_2d, xedges, yedges = np.histogram2d(
        x_data, y_data, range=[[xmin, xmax], [ymin, ymax]], bins=binsxy)

    hist = gaussian_filter(hist_2d, std_dev, mode='nearest')  # 'constant')

    return hist, xedges, yedges


def detect_peaks(image):
    """
    Takes an image and detect the peaks using the local maximum filter.
    Returns a boolean mask of the peaks (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """

    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2, 2)

    # apply the local maximum filter; all pixel of maximal value
    # in their neighborhood are set to 1
    local_max = maximum_filter(image, footprint=neighborhood) == image
    # local_max is a mask that contains the peaks we are
    # looking for, but also the background.
    # In order to isolate the peaks we must remove the background from the
    # mask.

    # we create the mask of the background
    background = (image == 0)

    # a little technicality: we must erode the background in order to
    # successfully subtract it form local_max, otherwise a line will
    # appear along the background border (artifact of the local maximum filter)
    eroded_background = binary_erosion(
        background, structure=neighborhood, border_value=1)

    # we obtain the final mask, containing only peaks,
    # by removing the background from the local_max mask
    detected_peaks = local_max - eroded_background

    return detected_peaks


def get_peaks_coords(peaks, ra, dec, histo_edges):
    '''
    Transform peak bin coordinates into actual coordinates.
    '''

    peaks_coords = []
    # for each smoothed chart.
    for i, h_p in enumerate(peaks):
        xedges, yedges = histo_edges[i]
        max_dens_coords = []
        # For each x column in the histogram.
        for x_p_indx, x_col in enumerate(h_p):
            # For each y column in the histogram.
            for y_p_indx, y_p in enumerate(x_col):
                if y_p:
                    x_cent_bin, y_cent_bin = x_p_indx, y_p_indx
                    # x,y coords of the center of the above bin.
                    x_cent_pix = np.average(xedges[x_cent_bin:x_cent_bin + 2])
                    y_cent_pix = np.average(yedges[y_cent_bin:y_cent_bin + 2])
                    max_dens_coords.append([x_cent_pix, y_cent_pix])

        peaks_coords.append(max_dens_coords)

    return peaks_coords


def get_kde_dens(ra, dec, peaks_coords):
    '''
    '''

    # Obtain Gaussian KDE in coordinates.
    values = np.vstack([ra, dec])
    kernel = stats.gaussian_kde(values)

    max_dens = []
    for p in peaks_coords:
        # Evaluate kernel in peak coordinates.
        positions = np.vstack([zip(*p)[0], zip(*p)[1]])
        k_pos = kernel(positions)
        # Store density values in those positions.
        max_dens.append(k_pos)

    return max_dens


def get_filter_dens(peaks_coords, max_dens, thresh=0.75):
    '''
    '''

    max_dens_coords, max_dens_filt = [], []
    for p, md in zip(peaks_coords, max_dens):
        max_d = max(md)
        coords_f, dens_f = [], []
        for i, dens in enumerate(md):

            # Only store those that are larger than a certain threshold.
            if dens > (max_d * thresh):
                dens_f.append(dens)
                coords_f.append(p[i])

        max_dens_coords.append(coords_f)
        max_dens_filt.append(dens_f)

    return max_dens_coords, max_dens_filt


def make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
              max_dens_coords, max_dens_filt, thresh):
    '''
    '''

    fig = plt.figure(figsize=(15, 50))
    gs = gridspec.GridSpec(16, 3)

    i, j = 0, 0
    for h, p, m, mc, mf in zip(clust_histo, peaks_coords, max_dens,
                               max_dens_coords, max_dens_filt):
        print 'Plotting: ', i, i + 1, i + 2

        ax1 = plt.subplot(gs[i])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.imshow(h.transpose(), origin='lower', interpolation='none')
        # ax.imshow(h.transpose(), origin='lower')
        # Text box.
        ob = offsetbox.AnchoredText(r'$\sigma={}$'.format(std_dev_lst[j]),
                                    loc=2, prop=dict(size=8))
        ob.patch.set(alpha=0.65)
        ax1.add_artist(ob)
        ax1.invert_xaxis()
        j += 1

        # ax2 = plt.subplot(gs[i + 1])
        # cm = plt.cm.gist_earth_r
        # plt.imshow(p.transpose(), origin='lower', cmap=cm)
        # ax2.invert_xaxis()

        siz = 5 + 95 * (np.array(m) / max(m)) ** 2
        ax2 = plt.subplot(gs[i + 1])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
        plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
        x_c, y_c = zip(*p)[0], zip(*p)[1]
        plt.xlim(min(ra), max(ra))
        plt.ylim(min(dec), max(dec))
        cm = plt.cm.get_cmap('RdYlBu_r')
        plt.scatter(x_c, y_c, c=m, cmap=cm, s=siz, lw=0.25)
        ax2.invert_xaxis()
        ax2.set_aspect('equal')

        # v_min_mp, v_max_mp = min(m), max(m)
        siz = 5 + 95 * (np.array(mf) / max(mf)) ** 6
        ax3 = plt.subplot(gs[i + 2])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
        plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
        x_c, y_c = zip(*mc)[0], zip(*mc)[1]
        plt.xlim(min(ra), max(ra))
        plt.ylim(min(dec), max(dec))
        cm = plt.cm.get_cmap('RdYlBu_r')
        plt.scatter(x_c, y_c, c=mf, cmap=cm, s=siz, lw=0.25)
        # vmin=v_min_mp, vmax=v_max_mp)
        # Text box.
        ob = offsetbox.AnchoredText('thresh={}'.format(thresh), loc=2,
                                    prop=dict(size=8))
        ob.patch.set(alpha=0.5)
        ax3.add_artist(ob)
        # Invert axis and set aspect.
        ax3.invert_xaxis()
        ax3.set_aspect('equal')

        # Increase counter.
        i += 3

    # Output png file.
    fig.tight_layout()
    plt.savefig('/home/gabriel/Descargas/peak_detect.png', dpi=300)


# Get coordinates.
ra, dec = get_coords()

# Generate 2D histogram of the coordinates, using several different numbers
# of bins and standard deviating values for the Gaussian smoothing.
clust_histo, histo_edges = [], []
std_dev_lst = [2, 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 4.25, 4.5,
               4.75, 5., 5.25, 5.5, 5.75]
for N_bins in [100]:
    for std_dev in std_dev_lst:
        h, xedges, yedges = get_2d_histo(ra, dec, N_bins, std_dev)
        clust_histo.append(h)
        histo_edges.append([xedges, yedges])

# Detect peaks.
peaks = []
for h in clust_histo:
    peaks.append(detect_peaks(h))

# Transform peaks bin coordinates into actual coordinates.
peaks_coords = get_peaks_coords(peaks, ra, dec, histo_edges)

# Obtain KDE values for each peak in each smoothed cluster.
max_dens = get_kde_dens(ra, dec, peaks_coords)

thresh = 0.9
# Filter found max density coordinates below a certain threshold.
max_dens_coords, max_dens_filt = get_filter_dens(peaks_coords, max_dens,
                                                 thresh)

# Plot results.
make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
          max_dens_coords, max_dens_filt, thresh)

	import numpy as np
	from scipy.ndimage.filters import maximum_filter
	from scipy.ndimage.filters import gaussian_filter
	from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
	from scipy import stats
	import matplotlib.pyplot as plt
	import matplotlib.gridspec as gridspec
	import matplotlib.offsetbox as offsetbox
	# from mpl_toolkits.axes_grid1 import make_axes_locatable


	def get_coords():
	'''
	Return coordinates from file.
	'''
	# Each sub-list in file_data is a row of the file.
	# file_data = np.loadtxt("TR24_clean.dat")
	file_data = np.loadtxt("CLUSTER.DAT")

	# Extract coordinates and zip them into two lists.
	ra, dec = zip(file_data)[1], zip(file_data)[2]

	# ra_rang, dec_rang = max(ra) - min(ra), max(dec) - min(dec)
	# print 'RA range: ', ra_rang
	# print 'DEC range: ', dec_rang

	return ra, dec


	def get_2d_histo(x_data, y_data, N_bins, std_dev):
	'''
	Return 2D histogram for the positional data.
	'''
	xmin, xmax = min(x_data), max(x_data)
	ymin, ymax = min(y_data), max(y_data)

	# Calculate the number of bins used.
	x_rang, y_rang = (xmax - xmin), (ymax - ymin)

	# Bin width to create the 2D histogram.
	bin_width = min(x_rang, y_rang) / N_bins

	# Number of bins in x,y given the bin width.
	binsxy = [int(x_rang / bin_width), int(y_rang / bin_width)]

	# hist_2d is the 2D histogram, *edges store the edges of the bins.
	hist_2d, xedges, yedges = np.histogram2d(
	x_data, y_data, range=[[xmin, xmax], [ymin, ymax]], bins=binsxy)

	hist = gaussian_filter(hist_2d, std_dev, mode='nearest') # 'constant')

	return hist, xedges, yedges


	def detect_peaks(image):
	"""
	Takes an image and detect the peaks using the local maximum filter.
	Returns a boolean mask of the peaks (i.e. 1 when
	the pixel's value is the neighborhood maximum, 0 otherwise)
	"""

	# define an 8-connected neighborhood
	neighborhood = generate_binary_structure(2, 2)

	# apply the local maximum filter; all pixel of maximal value
	# in their neighborhood are set to 1
	local_max = maximum_filter(image, footprint=neighborhood) == image
	# local_max is a mask that contains the peaks we are
	# looking for, but also the background.
	# In order to isolate the peaks we must remove the background from the
	# mask.

	# we create the mask of the background
	background = (image == 0)

	# a little technicality: we must erode the background in order to
	# successfully subtract it form local_max, otherwise a line will
	# appear along the background border (artifact of the local maximum filter)
	eroded_background = binary_erosion(
	background, structure=neighborhood, border_value=1)

	# we obtain the final mask, containing only peaks,
	# by removing the background from the local_max mask
	detected_peaks = local_max - eroded_background

	return detected_peaks


	def get_peaks_coords(peaks, ra, dec, histo_edges):
	'''
	Transform peak bin coordinates into actual coordinates.
	'''

	peaks_coords = []
	# for each smoothed chart.
	for i, h_p in enumerate(peaks):
	xedges, yedges = histo_edges[i]
	max_dens_coords = []
	# For each x column in the histogram.
	for x_p_indx, x_col in enumerate(h_p):
	# For each y column in the histogram.
	for y_p_indx, y_p in enumerate(x_col):
	if y_p:
	x_cent_bin, y_cent_bin = x_p_indx, y_p_indx
	# x,y coords of the center of the above bin.
	x_cent_pix = np.average(xedges[x_cent_bin:x_cent_bin + 2])
	y_cent_pix = np.average(yedges[y_cent_bin:y_cent_bin + 2])
	max_dens_coords.append([x_cent_pix, y_cent_pix])

	peaks_coords.append(max_dens_coords)

	return peaks_coords


	def get_kde_dens(ra, dec, peaks_coords):
	'''
	'''

	# Obtain Gaussian KDE in coordinates.
	values = np.vstack([ra, dec])
	kernel = stats.gaussian_kde(values)

	max_dens = []
	for p in peaks_coords:
	# Evaluate kernel in peak coordinates.
	positions = np.vstack([zip(p)[0], zip(p)[1]])
	k_pos = kernel(positions)
	# Store density values in those positions.
	max_dens.append(k_pos)

	return max_dens


	def get_filter_dens(peaks_coords, max_dens, thresh=0.75):
	'''
	'''

	max_dens_coords, max_dens_filt = [], []
	for p, md in zip(peaks_coords, max_dens):
	max_d = max(md)
	coords_f, dens_f = [], []
	for i, dens in enumerate(md):

	# Only store those that are larger than a certain threshold.
	if dens > (max_d * thresh):
	dens_f.append(dens)
	coords_f.append(p[i])

	max_dens_coords.append(coords_f)
	max_dens_filt.append(dens_f)

	return max_dens_coords, max_dens_filt


	def make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
	max_dens_coords, max_dens_filt, thresh):
	'''
	'''

	fig = plt.figure(figsize=(15, 50))
	gs = gridspec.GridSpec(16, 3)

	i, j = 0, 0
	for h, p, m, mc, mf in zip(clust_histo, peaks_coords, max_dens,
	max_dens_coords, max_dens_filt):
	print 'Plotting: ', i, i + 1, i + 2

	ax1 = plt.subplot(gs[i])
	plt.tick_params(axis='both', which='major', labelsize=8)
	plt.imshow(h.transpose(), origin='lower', interpolation='none')
	# ax.imshow(h.transpose(), origin='lower')
	# Text box.
	ob = offsetbox.AnchoredText(r'$\sigma={}$'.format(std_dev_lst[j]),
	loc=2, prop=dict(size=8))
	ob.patch.set(alpha=0.65)
	ax1.add_artist(ob)
	ax1.invert_xaxis()
	j += 1

	# ax2 = plt.subplot(gs[i + 1])
	# cm = plt.cm.gist_earth_r
	# plt.imshow(p.transpose(), origin='lower', cmap=cm)
	# ax2.invert_xaxis()

	siz = 5 + 95 * (np.array(m) / max(m)) ** 2
	ax2 = plt.subplot(gs[i + 1])
	plt.tick_params(axis='both', which='major', labelsize=8)
	plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
	plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
	x_c, y_c = zip(p)[0], zip(p)[1]
	plt.xlim(min(ra), max(ra))
	plt.ylim(min(dec), max(dec))
	cm = plt.cm.get_cmap('RdYlBu_r')
	plt.scatter(x_c, y_c, c=m, cmap=cm, s=siz, lw=0.25)
	ax2.invert_xaxis()
	ax2.set_aspect('equal')

	# v_min_mp, v_max_mp = min(m), max(m)
	siz = 5 + 95 * (np.array(mf) / max(mf)) ** 6
	ax3 = plt.subplot(gs[i + 2])
	plt.tick_params(axis='both', which='major', labelsize=8)
	plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
	plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
	x_c, y_c = zip(mc)[0], zip(mc)[1]
	plt.xlim(min(ra), max(ra))
	plt.ylim(min(dec), max(dec))
	cm = plt.cm.get_cmap('RdYlBu_r')
	plt.scatter(x_c, y_c, c=mf, cmap=cm, s=siz, lw=0.25)
	# vmin=v_min_mp, vmax=v_max_mp)
	# Text box.
	ob = offsetbox.AnchoredText('thresh={}'.format(thresh), loc=2,
	prop=dict(size=8))
	ob.patch.set(alpha=0.5)
	ax3.add_artist(ob)
	# Invert axis and set aspect.
	ax3.invert_xaxis()
	ax3.set_aspect('equal')

	# Increase counter.
	i += 3

	# Output png file.
	fig.tight_layout()
	plt.savefig('/home/gabriel/Descargas/peak_detect.png', dpi=300)


	# Get coordinates.
	ra, dec = get_coords()

	# Generate 2D histogram of the coordinates, using several different numbers
	# of bins and standard deviating values for the Gaussian smoothing.
	clust_histo, histo_edges = [], []
	std_dev_lst = [2, 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 4.25, 4.5,
	4.75, 5., 5.25, 5.5, 5.75]
	for N_bins in [100]:
	for std_dev in std_dev_lst:
	h, xedges, yedges = get_2d_histo(ra, dec, N_bins, std_dev)
	clust_histo.append(h)
	histo_edges.append([xedges, yedges])

	# Detect peaks.
	peaks = []
	for h in clust_histo:
	peaks.append(detect_peaks(h))

	# Transform peaks bin coordinates into actual coordinates.
	peaks_coords = get_peaks_coords(peaks, ra, dec, histo_edges)

	# Obtain KDE values for each peak in each smoothed cluster.
	max_dens = get_kde_dens(ra, dec, peaks_coords)

	thresh = 0.9
	# Filter found max density coordinates below a certain threshold.
	max_dens_coords, max_dens_filt = get_filter_dens(peaks_coords, max_dens,
	thresh)

	# Plot results.
	make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
	max_dens_coords, max_dens_filt, thresh)