Simple Point Source Example

This tutorial shows you hwo to work with pylira using a simple simulated point source and Gaussian shaped point spread function.

We first define the required imports:

%matplotlib inline
import copy
import numpy as np
import matplotlib.pyplot as plt
from regions import CirclePixelRegion, PixCoord
from multiprocessing import Pool

from pylira.data import point_source_gauss_psf, gauss_and_point_sources_gauss_psf
from pylira.utils.plot import plot_example_dataset, plot_hyperpriors_lira
from pylira import LIRADeconvolver, LIRADeconvolverResult, LIRASignificanceEstimator

random_state = np.random.RandomState(148)

Dataset

Then we can use the method point_source_gauss_psf to get a pre-defined test dataset:

data = point_source_gauss_psf(random_state=random_state)
print(data.keys())

dict_keys(['counts', 'psf', 'exposure', 'background', 'flux'])

The data variable is a dict containing the counts, psf, exposure, background and ground truth for flux. We can illustrate the data using:

plot_example_dataset(data)

Configure and Run Deconvolution

Next we define the LIRADeconvolver class, which holds the configuration the algorithm will be run with:

alpha_init = 0.1 * np.ones(np.log2(data["counts"].shape[0]).astype(int))

dec = LIRADeconvolver(
    alpha_init=alpha_init,
    n_iter_max=5_000,
    n_burn_in=200,
    fit_background_scale=False,
    random_state=random_state
)

We can print the instance to see the full configuration:

print(dec)

LIRADeconvolver
---------------

    alpha_init           : [0.1, 0.1, 0.1, 0.1, 0.1]
    n_iter_max           : 5000
    n_burn_in            : 200
    fit_background_scale : False
    save_thin            : True
    ms_ttlcnt_pr         : 1
    ms_ttlcnt_exp        : 0.05
    ms_al_kap1           : 0.0
    ms_al_kap2           : 1000.0
    ms_al_kap3           : 3.0

Now we can also visualize the hyperprior distribution of the alpha parameter:

plot_hyperpriors_lira(
    ncols=1,
    figsize=(8, 5),
    alphas=np.linspace(0, 0.3, 100),
    ms_al_kap1=0,
    ms_al_kap2=1000,
    ms_al_kap3=3
);

Now we have to define the initial guess for the flux and add it to the data dictionary. For this we use the reserved "flux_init"key:

data["flux_init"] = random_state.gamma(10, size=data["counts"].shape)

Finally we can run the LIRA algorithm using .run():

%%time
result = dec.run(data)

CPU times: user 1min 50s, sys: 469 ms, total: 1min 50s
Wall time: 1min 50s

Results and Diagnosis

To check the validity of the result we can plot the posterior mean:

result.plot_posterior_mean(from_image_trace=True)

As well as the parameter traces:

result.plot_parameter_traces()

/home/docs/checkouts/readthedocs.org/user_builds/pylira/envs/latest/lib/python3.7/site-packages/numpy/core/_methods.py:230: RuntimeWarning: invalid value encountered in subtract
  x = asanyarray(arr - arrmean)

We can also plot the parameter distributions:

result.plot_parameter_distributions()

And pixel traces in a given region, to check for correlations with neighbouring pixels:

result.plot_pixel_traces_region(center_pix=(16, 16), radius_pix=2);

Significance Estimation

The signifiance estimation follows the method proposed by Stein et al. 2015.

Null Hypothesis Simulations

First we simulate (“fake”) observations for the null hypothesis: in this case we assume the null hypothesis includes the background only:

# number of simulations
n_iter = 10

# reduce number of iterations for LIRA
dec.n_iter_max = 1_000

datasets = []
data_null = copy.deepcopy(data)

for idx in range(n_iter):
    data_null["counts"] = random_state.poisson(data["background"])
    data_null["flux_init"] = random_state.gamma(10, size=data["counts"].shape)
    datasets.append(data_null)

%%time
with Pool(processes=8) as pool:
    results = pool.map(dec.run, datasets)

CPU times: user 259 ms, sys: 137 ms, total: 396 ms
Wall time: 2min 45s

Define Regions of Interest

# Region of interest for the source
region_src = CirclePixelRegion(
    center=PixCoord(16, 16),
    radius=2
)

# Some control region
region_bkg = CirclePixelRegion(
    center=PixCoord(23, 8),
    radius=5
)

labels = np.zeros((32, 32))

for idx, region in enumerate([region_src, region_bkg]):
    mask = region.to_mask()
    labels += (idx + 1) * mask.to_image(shape=labels.shape)

plt.imshow(labels)
plt.colorbar();

Estimate and Visualize Distributions of \(\xi\)

First we instantiate the LIRASignificanceEstimator with the result obtained before, the results from the null hypothesis simulations and the label image:

est = LIRASignificanceEstimator(
    result_observed_im=result,
    result_replicates=results,
    labels_im=labels
)

Now we first estimate the p-values for the regions defined by the labels.

result_p_values = est.estimate_p_values(
    data=data, gamma=0.005
)

For the we consider the parameter \(\xi\), defined as:

\(\xi = \tau_1 / (\tau_1 + \tau_0)\)

So the proportion of the total image intensity that is due to the added component \(\tau_1\). The distribution of \(\xi\) can be plotted the different regions. First the region of interest in the center:

ax = est.plot_xi_dist(
    xi_obs=result_p_values[2],
    xi_repl=result_p_values[1],
    region_id="1.0"
);
ax.set_ylim(0, 0.2);
ax.set_xlim(-8.3, 0.2);

In this case all the probability mass is in the region of the bright point source and the mean of \(\xi\) approaches 1.

We also illustrate the distribution of \(\xi\) for second region, where we only expect background:

ax = est.plot_xi_dist(
    xi_obs=result_p_values[2],
    xi_repl=result_p_values[1],
    region_id="2.0"
);
ax.set_ylim(0, 0.2);
ax.set_xlim(-8.3, 0.2);

In this case the distribution is fully compatible with the distribution of the null hypothesis.

Notebook file for download: point_source.ipynb