<commit_before>#!/usr/bin/env python
# encoding: utf-8
"""
deconvolve.py

GALFIT "deconvolution" as in Szomoru et al.
https://iopscience.iop.org/article/10.1088/2041-8205/714/2/L244

"""

import numpy as np

import astropy.io.fits as pyfits
import astropy.stats
import astropy.table


def test():

    im = pyfits.open('/tmp/gf_out.fits')
    model_hdu = im[2]
    residuals = im[3].data

    _x = galfit_deconvolve(model_hdu, residuals, rms=None)

    model_hdu = _x['model']
    residuals = _x['resid'].data

    imshow_args = {'vmin': -5, 'vmax': 5, 'cmap': 'viridis'}

    _ = galfit_deconvolve(model_hdu, residuals, rms=_x['rms'], imshow_args=imshow_args)


def galfit_deconvolve(model_hdu, residuals, rms=None, mask=None, oversample=8, comp='1', xprof=np.append(0, np.arange(2, 20)), y_func=np.mean, cumul_values=[0.5, 0.8], make_plot=True, imshow_args={'vmin': -5, 'vmax': 5, 'cmap': 'viridis'}, plt_kwargs={'linestyle': 'steps-mid', 'color': 'r', 'alpha': 0.8}, npix=20, psf_offset=[1, 1]):
    """
    Deconvolve an image using a galfit Sersic model

    `model_hdu` is the HDU output from GALFIT containing the model and model
    parameters (e.g, out.fits[2]).

    `residuals` are the fit residuals (out.fits[3].data)

    `psf_offset` is the offset of the PSF relative to the central pixel of
    the PSF array.

    Returns:

        R = radii evaluated at pixel positions for the ellipse parameters
        phi = position angle around the ellipse, radians
        prof = 2D array of unconvolved Sersic profile
        params = Galfit params
        tab = table of the 1D profile
        fig = optional figure if `make_plot=True`

    """
    import matplotlib.pyplot as plt

    #model_hdu = _x['model']
    #residuals = _x['resid']

    _h = model_hdu.header
    shape = residuals.shape

    if mask is None:
        mask = np.isfinite(residuals)

    # Ellipse parameters and sersic profile
    R, phi, prof, params = sersic_profile(shape, oversample=oversample,
                                          gf_header=_h, comp=comp,
                                          psf_offset=psf_offset)

    so = np.argsort(R[mask].flatten())
    Rso = R[mask].flatten()[so]

    data = residuals + model_hdu.data

    if '1_SKY' in _h:
        sky = gf_header_key(_h, '1_SKY')
    else:
        sky = 0.

    ##########
    # 1D averaged profiles

    integ = False

    # Original data
    _, ydata, _, _ = running_median(Rso,
                (data)[mask].flatten()[so].astype(np.float)-sky,
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)

    # Sersic profile
    _, yprof, _, _ = running_median(Rso,
                (prof)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)

    # Sersic + residuals
    xm, ydeconv, ys, yn = running_median(Rso,
                (prof + residuals)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)

    # Sum for total normalizatioin
    _, ydeconv_sum, _, _ = running_median(Rso,
                (prof + residuals)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=np.sum, integrate=False)

    # Variance
    if rms is not None:
        xmx, yv, yvs, yn = running_median(Rso,
                rms[mask].flatten()[so].astype(np.float)**2,
                bins=xprof, use_median=False, y_func=np.sum)

        yrms = np.sqrt(yv)/yn
        im_norm = rms
    else:
        yrms = ys
        im_norm = 1

    dx = np.diff(xprof)
    xpix = xprof[1:]-dx/2.

    if False:
        # Scale
        yscl = 1.

        #yscl = prof[msk].sum()/np.trapz(2*np.pi*yprof*xpix, xpix)
        yprof *= yscl
        ydeconv *= yscl
        yrms *= yscl

    # Interpolate Radii from the cumulative flux distribution
    cumul_flux = np.cumsum((prof + residuals)[mask].flatten()[so])

    #total = cumflux.max()
    total = ydeconv_sum.sum()

    Rcumul = np.interp(cumul_values, cumul_flux/total, Rso)

    tab = astropy.table.Table()
    tab['xpix'] = xpix
    tab['yprofile'] = yprof
    tab['ydeconv'] = ydeconv
    tab['yrms'] = yrms
    tab['ydata'] = ydata

    tab.meta['total_flux'] = total
    for k in params:
        tab.meta[k] = params[k]

    for i, r in enumerate(cumul_values):
        tab.meta['R{0}'.format(int(r*100))] = Rcumul[i]

    params['Rcumul'] = Rcumul

    if make_plot:

        fig = plt.figure(figsize=(12, 3))
        ax = fig.add_subplot(141)
        ax2 = fig.add_subplot(142)
        ax3 = fig.add_subplot(143)<commit_msg>add model to galfit deconvolved<commit_after>#!/usr/bin/env python
# encoding: utf-8
"""
deconvolve.py

GALFIT "deconvolution" as in Szomoru et al.
https://iopscience.iop.org/article/10.1088/2041-8205/714/2/L244

"""

import numpy as np

import astropy.io.fits as pyfits
import astropy.stats
import astropy.table


def test():

    im = pyfits.open('/tmp/gf_out.fits')
    model_hdu = im[2]
    residuals = im[3].data

    _x = galfit_deconvolve(model_hdu, residuals, rms=None)

    model_hdu = _x['model']
    residuals = _x['resid'].data

    imshow_args = {'vmin': -5, 'vmax': 5, 'cmap': 'viridis'}

    _ = galfit_deconvolve(model_hdu, residuals, rms=_x['rms'], imshow_args=imshow_args)


def galfit_deconvolve(model_hdu, residuals, rms=None, mask=None, oversample=8, comp='1', xprof=np.append(0, np.arange(2, 20)), y_func=np.mean, cumul_values=[0.5, 0.8], make_plot=True, imshow_args={'vmin': -5, 'vmax': 5, 'cmap': 'viridis'}, plt_kwargs={'linestyle': 'steps-mid', 'color': 'r', 'alpha': 0.8}, npix=20, psf_offset=[1, 1]):
    """
    Deconvolve an image using a galfit Sersic model

    `model_hdu` is the HDU output from GALFIT containing the model and model
    parameters (e.g, out.fits[2]).

    `residuals` are the fit residuals (out.fits[3].data)

    `psf_offset` is the offset of the PSF relative to the central pixel of
    the PSF array.

    Returns:

        R = radii evaluated at pixel positions for the ellipse parameters
        phi = position angle around the ellipse, radians
        prof = 2D array of unconvolved Sersic profile
        params = Galfit params
        tab = table of the 1D profile
        fig = optional figure if `make_plot=True`

    """
    import matplotlib.pyplot as plt

    #model_hdu = _x['model']
    #residuals = _x['resid']

    _h = model_hdu.header
    shape = residuals.shape

    if mask is None:
        mask = np.isfinite(residuals)

    # Ellipse parameters and sersic profile
    R, phi, prof, params = sersic_profile(shape, oversample=oversample,
                                          gf_header=_h, comp=comp,
                                          psf_offset=psf_offset)

    so = np.argsort(R[mask].flatten())
    Rso = R[mask].flatten()[so]

    data = residuals + model_hdu.data

    if '1_SKY' in _h:
        sky = gf_header_key(_h, '1_SKY')
    else:
        sky = 0.

    ##########
    # 1D averaged profiles

    integ = False

    # Original data
    _, ydata, _, _ = running_median(Rso,
                (data)[mask].flatten()[so].astype(np.float)-sky,
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)
    
    # Convolved model
    _, ymodel, _, _ = running_median(Rso,
                (model_hdu.data)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)
    
    # Sersic profile
    _, yprof, _, _ = running_median(Rso,
                (prof)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)

    # Sersic + residuals
    xm, ydeconv, ys, yn = running_median(Rso,
                (prof + residuals)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=y_func, integrate=integ)

    # Sum for total normalizatioin
    _, ydeconv_sum, _, _ = running_median(Rso,
                (prof + residuals)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=np.sum, integrate=False)

    # Variance
    if rms is not None:
        xmx, yv, yvs, yn = running_median(Rso,
                rms[mask].flatten()[so].astype(np.float)**2,
                bins=xprof, use_median=False, y_func=np.sum)

        yrms = np.sqrt(yv)/yn
        im_norm = rms
        
        # weighted
        xmx, ynum_model, yvs, yn = running_median(Rso,
                (model_hdu.data/rms**2)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=np.sum)

        xmx, ynum, yvs, yn = running_median(Rso,
                (data/rms**2)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=np.sum)
                
        xmx, yden, yvs, yn = running_median(Rso,
                (1./rms**2)[mask].flatten()[so].astype(np.float),
                bins=xprof, use_median=False, y_func=np.sum)
        
        yweight = ynum/yden
        yweight_model = ynum_model/yden
        yweight_err = 1./np.sqrt(yden)
        
    else:
        yrms = ys
        im_norm = 1
        yweight = None
        
    dx = np.diff(xprof)
    xpix = xprof[1:]-dx/2.

    if False:
        # Scale
        yscl = 1.

        #yscl = prof[msk].sum()/np.trapz(2*np.pi*yprof*xpix, xpix)
        yprof *= yscl
        ydeconv *= yscl
        yrms *= yscl

    # Interpolate Radii from the cumulative flux distribution
    cumul_flux = np.cumsum((prof + residuals)[mask].flatten()[so])

    #total = cumflux.max()
    total = ydeconv_sum.sum()

    Rcumul = np.interp(cumul_values, cumul_flux/total, Rso)

    tab = astropy.table.Table()
    tab['xpix'] = xpix
    tab['yprofile'] = yprof
    tab['ydeconv'] = ydeconv
    tab['yrms'] = yrms
    tab['ydata'] = ydata
    tab['ymodel'] = ymodel
    
    if yweight is not None:
        tab['yweight'] = yweight
        tab['yweight_model'] = yweight_model
        tab['yweight_err'] = yweight_err
        
    tab.meta['total_flux'] = total
    for k in params:
        tab.meta[k] = params[k]

    for i, r in enumerate(cumul_values):
        tab.meta['R{0}'.format(int(r*100))] = Rcumul[i]

    params['Rcumul'] = Rcumul

    if make_plot:

        fig = plt.figure(figsize=(12, 3))
        ax = fig.add_subplot(141)
        ax2 = fig.add_subplot(142)
        ax3 = fig.add_subplot(143)
{'fraction_non_alphanumeric': 0.13466550825369245, 'fraction_numerical': 0.015155903079447823, 'mean_word_length': 2.4430043203722165, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5777858', 'n_tokens_mistral': 4016, 'n_tokens_neox': 3695, 'n_words': 785}