evaluation.py 12.6 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
from utils import find_nearest
LUSTIG Peter's avatar
LUSTIG Peter committed
30 31


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

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

plt.ion()


LUSTIG Peter's avatar
LUSTIG Peter committed
44 45 46 47
class PCEvaluation:
    def __init__(self, sources, fake_sources, shape, wcs, flux=None,
                 mapbins=19, threshold_bins=5, threshold_range=(3, 5)):

LUSTIG Peter's avatar
LUSTIG Peter committed
48
        print('sorting values')
LUSTIG Peter's avatar
LUSTIG Peter committed
49 50
        idxsort = np.argsort(flux.to_value(u.mJy))
        self.flux = flux[idxsort]
LUSTIG Peter's avatar
LUSTIG Peter committed
51 52
        self.sources = [sources[i] for i in idxsort]
        self.fake_sources = [fake_sources[i] for i in idxsort]
LUSTIG Peter's avatar
LUSTIG Peter committed
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

        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)
LUSTIG Peter's avatar
LUSTIG Peter committed
74 75
        #print(self.wcs_3D)
        #sys.exit()
LUSTIG Peter's avatar
LUSTIG Peter committed
76

LUSTIG Peter's avatar
LUSTIG Peter committed
77
        print('wcs created')
LUSTIG Peter's avatar
LUSTIG Peter committed
78 79 80 81 82
        # 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)
LUSTIG Peter's avatar
LUSTIG Peter committed
83 84
        self.completness, self.purity, self.hitmap = self.GetCP()
        print(self.completness.shape)
LUSTIG Peter's avatar
LUSTIG Peter committed
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104

    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)):
LUSTIG Peter's avatar
LUSTIG Peter committed
105
                tmpres = self.completness_purity(sources[i], fake_sources[i],
LUSTIG Peter's avatar
LUSTIG Peter committed
106 107 108 109
                                                 wcs, shape)
                comp.append(tmpres[0])
                pur.append(tmpres[1])
                hitm.append(tmpres[2])
LUSTIG Peter's avatar
LUSTIG Peter committed
110
            return np.array(comp), np.array(pur), np.array(hitm)
LUSTIG Peter's avatar
LUSTIG Peter committed
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

    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
142 143


LUSTIG Peter's avatar
LUSTIG Peter committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    def GetBinResults(self, ix, iy):
        #print(self.completness)
        return (self.completness[:, iy, ix, :],
                self.purity[:, iy, ix, :],
                self.hitmap[:, iy, ix])


    def PlotBin(self, data, title='', flux=np.array([]), thresh=[],
                       nfluxlabels=None, nthreshlabels=None, **kwargs):
        tickfs = 20
        labelfs = 25

        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]))

        plt.figure()
        plt.title(title, fontsize=30)
        plt.xlabel('Detection Threshold [SNR]', fontsize=labelfs)
        plt.ylabel('Flux [mJy]', fontsize=labelfs)
        plt.xticks(tlblpos, tlbl, fontsize=tickfs)
        plt.yticks(flblpos, flbl, fontsize=tickfs)
        plt.imshow(data, origin='lower', **kwargs)
        cbar = plt.colorbar()
        cbar.ax.tick_params(labelsize=tickfs)


190 191 192 193
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
194 195
sh = data.data.shape
wcs = data.wcs
196

LUSTIG Peter's avatar
LUSTIG Peter committed
197 198 199
hdul = fits.HDUList(fits.open('/home/peter/Dokumente/Uni/Paris/Stage/'
                              'FirstSteps/Completness/'
                              'combined_tables_long.fits'))
200 201 202
nfluxes = hdul[0].header['NFLUXES']
print('{} different fluxes found'.format(nfluxes))

LUSTIG Peter's avatar
LUSTIG Peter committed
203
indata = []
204

LUSTIG Peter's avatar
LUSTIG Peter committed
205 206 207 208
FLUX = []
SOURCE = []
FSOURCE = []

LUSTIG Peter's avatar
LUSTIG Peter committed
209
# for isimu in range(nfluxes):
LUSTIG Peter's avatar
LUSTIG Peter committed
210
for isimu in range(6):
LUSTIG Peter's avatar
LUSTIG Peter committed
211 212 213 214 215 216 217 218
    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])]))

xx = PCEvaluation(SOURCE, FSOURCE, sh, wcs, u.Quantity(FLUX))
LUSTIG Peter's avatar
LUSTIG Peter committed
219 220
print('done')
# %% testfunctions
LUSTIG Peter's avatar
LUSTIG Peter committed
221

LUSTIG Peter's avatar
LUSTIG Peter committed
222 223 224
# p = Pool(2)
# res = xx.GetCP()
dd = xx.GetBinResults(9, 9)
225

LUSTIG Peter's avatar
LUSTIG Peter committed
226 227 228 229 230
dd = xx.completness
print('comp shape', dd.shape)
xx.PlotBin(xx.completness[:, 9, 9, :])
plt.show(block=True)
sys.exit()
LUSTIG Peter's avatar
LUSTIG Peter committed
231
# %%
LUSTIG Peter's avatar
LUSTIG Peter committed
232
# midbin = int(bins / 2)
233

LUSTIG Peter's avatar
LUSTIG Peter committed
234 235 236

def PlotEvaluation(data, title='', flux=np.array([]), thresh=[],
                   nfluxlabels=None, nthreshlabels=None, **kwargs):
237 238
    tickfs = 20
    labelfs = 25
LUSTIG Peter's avatar
LUSTIG Peter committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262

    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]))

263 264 265 266
    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
267 268
    plt.xticks(tlblpos, tlbl, fontsize=tickfs)
    plt.yticks(flblpos, flbl, fontsize=tickfs)
269 270 271 272 273 274
    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
275 276 277
               flux=np.array(FLUX), thresh=threshold, cmap='bone',
               nfluxlabels=10, nthreshlabels=5, aspect='auto')

LUSTIG Peter's avatar
LUSTIG Peter committed
278 279 280 281 282
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
283 284 285 286 287

# %% 2D plot


def PlotFixedThreshold(thresholds, bin, completness, allthresholds, flux,
LUSTIG Peter's avatar
LUSTIG Peter committed
288
                       nfluxlabels=None, hlines=None):
LUSTIG Peter's avatar
LUSTIG Peter committed
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

    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)

315
# cmap bone hot
LUSTIG Peter's avatar
LUSTIG Peter committed
316 317 318 319 320 321
# 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
322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370

# %% 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)