"""Calculate PDF from set of maps and save output to HDF"""
import numpy as np
import pandas as pd
from astropy.io import fits


def cal_pdf(map_file, mask_file, bins=60):
    data = fits.getdata(map_file)[fits.getdata(mask_file).astype(bool)]
    pdf, bin_edges = np.histogram(data, bins=bins, density=True)
    bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])
    return pdf, bin_centers


if __name__ == '__main__':
    telescopes = ['hera37', 'hera61', 'hera91', 'hera127',
                  'hera169', 'hera217', 'hera271', 'hera331']
    prefixes = ['gauss_map', 'gauss_wm_imap']
    for tel in telescopes:
        for p in prefixes:
            maps_dir = '/data3/piyanat/runs/post_uvlt50_I_intnorm/maps/{:s}/'\
                .format(tel)
            outfile = '/data3/piyanat/runs/post_uvlt50_I_intnorm/pdf/' \
                      'mwa128_fhd_pdf_bw0.08MHz_windowing.h5'
            nbins = 60

            freqs, z, xi, ang = np.genfromtxt(
                '/data3/piyanat/runs/shared/'
                'interp_delta_21cm_f_z_xi_angpix.csv',
                delimiter=',', unpack=True
            )

            maps_files = [maps_dir + '{:s}_{:s}_{:.3f}MHz.fits'
                          .format(tel, p, f) for f in freqs]
            mask_files = [maps_dir + '{:s}_mask_map_{:.3f}MHz.fits'
                          .format(tel, f) for f in freqs]

            arr = np.empty((freqs.size, nbins, 2))

            for i in range(freqs.size):
                arr[i, :, 0], arr[i, :, 1] = cal_pdf(
                    maps_files[i], mask_files[i], bins=nbins
                )

            storage = pd.Panel(arr, items=['{:.3f}'.format(f) for f in freqs],
                               minor_axis=['pdf', 'bincen'])
            storage.to_hdf(outfile, 'pdf')
