evaluation.py 10.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
from __future__ import absolute_import, division, print_function

from pathlib import Path
import os
import numpy as np
import matplotlib.pyplot as plt

from multiprocessing import Pool, cpu_count
from functools import partial

from astropy import units as u
from astropy.io import ascii
from astropy.wcs import WCS
from astropy.utils.console import ProgressBar
from astropy.table import vstack

from scipy.optimize import curve_fit

from nikamap import NikaMap, Jackknife
from nikamap.utils import pos_uniform
from astropy.io import fits
from astropy.table import Table, MaskedColumn
import sys
from mpl_toolkits.axes_grid1 import make_axes_locatable
LUSTIG Peter's avatar
LUSTIG Peter committed
25 26 27
import dill as pickle
from matplotlib.ticker import FormatStrFormatter
from collections import OrderedDict
LUSTIG Peter's avatar
LUSTIG Peter committed
28
from utils import completness_purity_wcs, completness_worker, purity_worker
LUSTIG Peter's avatar
LUSTIG Peter committed
29 30


31 32 33 34 35 36 37 38 39 40 41 42

import os
os.getcwd()
'''
%load_ext autoreload
%autoreload 2
%matplotlib tk
'''

plt.ion()


LUSTIG Peter's avatar
LUSTIG Peter committed
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
class PCEvaluation:
    def __init__(self, sources, fake_sources, shape, wcs, flux=None,
                 mapbins=19, threshold_bins=5, threshold_range=(3, 5)):

        idxsort = np.argsort(flux.to_value(u.mJy))
        self.flux = flux[idxsort]
        self.sources = sources[idxsort]
        self.fake_sources = sources[idxsort]

        self.completness = None
        self.purity = None
        assert len(sources) == len(fake_sources), ("Number of results for "
                                                   "sources and fake "
                                                   "sources is not the same.")
        assert len(sources) == len(flux), ("Number of provided fluxes differs "
                                           "from number of simulation results")
        assert type(mapbins) is int, "number of bins must be an integer"
        self.bins = mapbins

        # threshold = np.linspace(threshold_range[0], threshold_range[1],
        #                         threshold_bins)
        threshold_edges = np.linspace(threshold_range[0], threshold_range[1],
                                      threshold_bins+1)

        self.shape_3D, self.wcs_3D = completness_purity_wcs(
                                        shape, wcs,
                                        bins=self.bins,
                                        threshold_range=threshold_range,
                                        threshold_bins=threshold_bins)

        # Testing the lower edges
        wcs_threshold = self.wcs_3D.sub([3])
        assert np.all(np.abs(wcs_threshold.all_pix2world(
                                np.arange(threshold_bins+1)-0.5, 0)
                             - threshold_edges) < 1e-15)

    def GetCP(self, sources=None, fake_sources=None, wcs=None, shape=None,
                pool=None):
        if sources is None:
            sources = self.sources
        if fake_sources is None:
            fake_sources = self.fake_sources
        if wcs is None:
            wcs = self.wcs_3D
        if shape is None:
            shape = self.shape_3D

        if pool is not None:
            f = partial(self.completness_purity, wcs=wcs, shape=shape)
            res = pool.starmap(f, (sources, fake_sources))
            res = list(zip(*res))
            return res[0], res[1], res[2]
        else:
            comp, pur, hitm = [], [], []
            for i in range(len(sources)):
                tmpres = self.completness_purity(sources, fake_sources,
                                                 wcs, shape)
                comp.append(tmpres[0])
                pur.append(tmpres[1])
                hitm.append(tmpres[2])
            return comp, pur, hitm

    def completness_purity(self, sources, fake_sources, wcs=None,
                           shape=None):
        """Compute completness map for a given flux"""

        min_threshold, max_threshold = wcs.sub([3]).all_pix2world(
                                                        [-0.5, shape[2]-1],
                                                        0)[0]

        # %load_ext snakeviz
        # %snakeviz the following line.... all is spend in the find_peaks /
        # fit_2d_gaussian
        # TODO: Change the find_peaks routine, or maybe just the
        # fit_2d_gaussian to be FAST ! (Maybe look into gcntrd.pro routine
        # or photutils.centroid.centroid_1dg maybe ?)

        completness, norm_comp = completness_worker(shape, wcs, sources,
                                                    fake_sources,
                                                    min_threshold,
                                                    max_threshold)

        purity, norm_pur = purity_worker(shape, wcs, sources, max_threshold)

        # norm can be 0, so to avoid warning on invalid values...
        with np.errstate(divide='ignore', invalid='ignore'):
            completness /= norm_comp[..., np.newaxis]
            purity /= norm_pur

        # TODO: One should probably return completness AND norm if one want to
        # combine several fluxes
        return completness, purity, norm_comp
LUSTIG Peter's avatar
LUSTIG Peter committed
135 136


137 138 139 140
plt.close('all')

DATA_DIR = "/home/peter/Dokumente/Uni/Paris/Stage/data/v_1"
data = NikaMap.read(Path(DATA_DIR) / '..' / 'map.fits')
LUSTIG Peter's avatar
LUSTIG Peter committed
141 142
sh = data.data.shape
wcs = data.wcs
143

LUSTIG Peter's avatar
LUSTIG Peter committed
144 145 146
hdul = fits.HDUList(fits.open('/home/peter/Dokumente/Uni/Paris/Stage/'
                              'FirstSteps/Completness/'
                              'combined_tables_long.fits'))
147 148 149
nfluxes = hdul[0].header['NFLUXES']
print('{} different fluxes found'.format(nfluxes))

LUSTIG Peter's avatar
LUSTIG Peter committed
150
indata = []
151

LUSTIG Peter's avatar
LUSTIG Peter committed
152 153 154 155
FLUX = []
SOURCE = []
FSOURCE = []

LUSTIG Peter's avatar
LUSTIG Peter committed
156
for isimu in range(nfluxes):
LUSTIG Peter's avatar
LUSTIG Peter committed
157 158 159 160 161 162 163 164 165 166
    FLUX.append(u.Quantity(hdul[0].header['flux{}'.format(isimu)]))

    SOURCE.append(Table.read(hdul['DETECTED_SOURCES{}'
                                  .format(FLUX[isimu])]))
    FSOURCE.append(Table.read(hdul['FAKE_SOURCES{}'
                                   .format(FLUX[isimu])]))

SOURCE = np.array(SOURCE)
FSOURCE = np.array(FSOURCE)
xx = PCEvaluation(SOURCE, FSOURCE, sh, wcs, u.Quantity(FLUX))
LUSTIG Peter's avatar
LUSTIG Peter committed
167

LUSTIG Peter's avatar
LUSTIG Peter committed
168
sys.exit()
169

LUSTIG Peter's avatar
LUSTIG Peter committed
170
# midbin = int(bins / 2)
171

LUSTIG Peter's avatar
LUSTIG Peter committed
172

173 174


LUSTIG Peter's avatar
LUSTIG Peter committed
175 176 177 178 179 180 181 182
def find_nearest(array, values):
    x, y = np.meshgrid(array, values)
    ev = np.abs(x - y)
    return np.argmin(ev, axis=1)


def PlotEvaluation(data, title='', flux=np.array([]), thresh=[],
                   nfluxlabels=None, nthreshlabels=None, **kwargs):
183 184
    tickfs = 20
    labelfs = 25
LUSTIG Peter's avatar
LUSTIG Peter committed
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210

    if nfluxlabels is not None:
        _label_flux = np.geomspace(flux[0], flux[-1], nfluxlabels)
        _f_idx = find_nearest(flux, _label_flux)
        flblpos, _flbl = _f_idx, flux[_f_idx]
    else:
        flblpos, _flbl = np.arange(len(flux)), flux

    if nthreshlabels is not None:
        _label_thresh = np.linspace(thresh[0], thresh[-1], nthreshlabels)
        _t_idx = find_nearest(thresh, _label_thresh)
        print(_t_idx)
        tlblpos, _tlbl = _t_idx, thresh[_t_idx]
    else:
        tlblpos, _tlbl = np.arange(len(thresh)), thresh

    flbl = []
    for i in range(len(_flbl)):
        flbl.append('{:.1f}'.format(_flbl[i]))

    tlbl = []
    for i in range(len(_tlbl)):
        tlbl.append('{:.1f}'.format(_tlbl[i]))
        # print(i, tlbl[i])

    # print(np.array([tlblpos, tlbl]).T)
211 212 213 214
    plt.figure()
    plt.title(title, fontsize=30)
    plt.xlabel('Detection Threshold [SNR]', fontsize=labelfs)
    plt.ylabel('Flux [mJy]', fontsize=labelfs)
LUSTIG Peter's avatar
LUSTIG Peter committed
215 216 217 218 219 220

    plt.xticks(tlblpos, tlbl, fontsize=tickfs)
    # ax = plt.gca()
    plt.yticks(flblpos, flbl, fontsize=tickfs)
    # ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))

221 222 223 224 225 226
    plt.imshow(data, origin='lower', **kwargs)
    cbar = plt.colorbar()
    cbar.ax.tick_params(labelsize=tickfs)


PlotEvaluation(COMPLETNESS[:, midbin, midbin, :], title='Completness',
LUSTIG Peter's avatar
LUSTIG Peter committed
227 228 229
               flux=np.array(FLUX), thresh=threshold, cmap='bone',
               nfluxlabels=10, nthreshlabels=5, aspect='auto')

LUSTIG Peter's avatar
LUSTIG Peter committed
230 231 232 233 234
PlotEvaluation(PURITY[:, midbin, midbin, :], title='Purity',
               flux=np.array(FLUX), thresh=threshold, cmap='bone',
               nfluxlabels=10, nthreshlabels=5, aspect='auto')
plt.show(block=True)

LUSTIG Peter's avatar
LUSTIG Peter committed
235 236 237 238 239

# %% 2D plot


def PlotFixedThreshold(thresholds, bin, completness, allthresholds, flux,
LUSTIG Peter's avatar
LUSTIG Peter committed
240
                       nfluxlabels=None, hlines=None):
LUSTIG Peter's avatar
LUSTIG Peter committed
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266

    linestyles = ['-', '--', '-.', ':']
    real_thresholds = find_nearest(allthresholds, thresholds)
    for i in range(len(real_thresholds)):
        _x = flux
        _y = completness[:, bin[0], bin[1], real_thresholds[i]]
        plt.plot(_x, _y, linestyle=linestyles[i],
                 label='{:.1f}'.format(allthresholds[real_thresholds[i]]))

    if hlines is not None:
        for i, val in enumerate(hlines):
            plt.axhline(val, color='r')
    plt.title('Fixed Threshold', fontsize=30, y=1.02)
    plt.xlabel('Source Flux [mJy]', fontsize=25)
    plt.ylabel('Completness', fontsize=25)
    plt.yticks(fontsize=20)
    plt.xticks(fontsize=20)
    plt.subplots_adjust(left=0.12)
    ax = plt.gca()
    ax.set_xscale("log", nonposx='clip')
    # legend = plt.legend(fontsize=25, title='SNR', loc='lower right')
    legend = plt.legend(fontsize=25, title='SNR', loc='upper left',
                        framealpha=1)
    plt.setp(legend.get_title(), fontsize=25)
    plt.show(block=True)

267
# cmap bone hot
LUSTIG Peter's avatar
LUSTIG Peter committed
268 269 270 271 272 273
# print(np.array(FLUX))


PlotFixedThreshold(np.array((3, 5, 7)), (midbin, midbin), COMPLETNESS,
                   threshold, np.array(FLUX), nfluxlabels=None,
                   hlines=[.5, .8, .9])
LUSTIG Peter's avatar
LUSTIG Peter committed
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322

# %% Plot map
# Plot_CompPur(completness, purity, threshold, nsim=None, savename=None,
#                  flux=None):

plotidx = np.array([24, 30, 40, 49])
# plotthreshidx = np.array([0, 4, 8])

plotcomp = COMPLETNESS[plotidx]
plotpur = PURITY[plotidx]
print(plotpur.shape)
# plotthresh = threshold[plotthreshidx]
plotflux = np.array(FLUX[plotidx])


def Plot_CompPur(completness, purity, threshold, nsim=None, savename=None,
                 flux=None):
    threshold_bins = completness.shape[-1]
    print('subplots:', threshold_bins)
    fig, axes = plt.subplots(nrows=2, ncols=threshold_bins, sharex=True,
                             sharey=True)
    print('axes shape', axes.shape)
    print('completness shape', completness.shape)
    for i in range(threshold_bins):
        print(i)
        axes[0, i].imshow(completness[:, :, i], vmin=0, vmax=1)
        im = axes[1, i].imshow(purity[:, :, i], vmin=0, vmax=1)
        axes[1, i].set_xlabel("thresh={:.2f}".format(threshold[i]))
        if i == (threshold_bins-1):
            # print('-----------')
            divider = make_axes_locatable(axes[1, i])
            cax = divider.append_axes('right', size='5%', pad=0.0)
            fig = plt.gcf()
            fig.colorbar(im, cax=cax, orientation='vertical')
    if nsim is not None:
        axes[0, 0].set_title("{} simulations".format(nsim))
    if flux is not None:
        axes[0, 1].set_title("{}".format(flux))
    axes[0, 0].set_ylabel("completness")
    axes[1, 0].set_ylabel("purity")
    if savename is not None:
        plt.savefig(savename)


for i in range(len(plotidx)):
    _comp = plotcomp[i]
    _pur = plotpur[i]
    Plot_CompPur(_comp, _pur, threshold, flux=plotflux[i])
plt.show(block=True)