Skip to content

Commit

Permalink
Merge pull request #483 from lmauviard/459-residuals-plotting-options
Browse files Browse the repository at this point in the history 
459 residuals plotting options
  • Loading branch information
@lmauviard
lmauviard authored Dec 16, 2024
2 parents 7b6df01 52a8486 commit f3a7e81
Show file tree
Hide file tree
Showing 12 changed files with 746 additions and 135 deletions.
307 changes: 264 additions & 43 deletions docs/source/Post-processing.ipynb
View file

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions examples/examples_fast/Outputs/ST_live_1000_eff_0.3_seed0_v2post_equal_weights.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 1 @@
0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 -0.554647675200076410E 05
2 changes: 2 additions & 0 deletions examples/examples_fast/Outputs/ST_live_1000_eff_0.3_seed0_v2summary.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 1,2 @@
0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 0.514165644197529786E-08 0.000000000000000000E 00 NaN 0.000000000000000000E 00 0.000000000000000000E 00 0.117746097693800428E-07 0.603268879392227826E-08 NaN 0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 -0.146886492504450813E 06 -0.554647675200076410E 05
0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 0.514165644197529786E-08 0.000000000000000000E 00 NaN 0.000000000000000000E 00 0.000000000000000000E 00 0.117746097693800428E-07 0.603268879392227826E-08 NaN 0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 0.140117499828338632E 01 0.116140460968017578E 02 0.908264352013451459E 00 0.720717012882232666E 00 -0.188007950782775879E-01 0.145097667350371728E 01 0.931132949382279529E 00 0.658978956937789917E 01 -0.554706595430130692E 05 -0.554647675200076410E 05
129 changes: 129 additions & 0 deletions xpsi/PostProcessing/_1d_residual.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 1,129 @@
from ._global_imports import *
from ._signalplot import SignalPlot
from scipy.stats import norm

class Residual1DPlot(SignalPlot):
""" Plot the count data, the posterior-expected count signal, and residuals.
The figure contain one panel :
* histogram of the standardised residuals between the
data and posterior-expected count signal over joint channel-phase
intervals.
The following example is (improved) from `Riley et al. 2019 <https://ui.adsabs.harvard.edu/abs/2019ApJ...887L..21R/abstract>`_:
.. image:: _static/_1d_residualplot.png
:param ... Todo
"""

__figtype__ = 'signalplot_1dresiduals'

# do not change at runtime (see base class comments):
__caching_targets__ = ['expected_counts']

__rows__ = 1
__columns__ = 1
__ax_rows__ = 1
__ax_columns__ = 1
__height_ratios__ = [1]
__width_ratios__ = [1] # second column for colorbars
__wspace__ = 0.025
__hspace__ = 0.125

@make_verbose('Instantiating a 1D residual plotter for residual distribution checking',
'Residual plotter instantiated')
def __init__(self,
nbins=20,
plot_fit=False,
parameters_vector=None,
**kwargs):
super(Residual1DPlot, self).__init__(**kwargs)

# Get the input parameters
self._nbins = int( nbins )
self._plot_fit = plot_fit
if not parameters_vector is None:
self.parameters_vector = parameters_vector

# Get on plotting
self._get_figure()
self._ax_1d_residuals = self._add_subplot(0,0)
self._ax_1d_residuals.set_xlabel(r'$(c_{ik}-d_{ik})/\sqrt{c_{ik}}$')
self._ax_1d_residuals.set_ylabel('Occurence')

plt.close()

@make_verbose('Residual1DPlot object iterating over samples',
'Residual1DPlot object finished iterating')
def execute(self, thetas, wrapper):
""" Loop over posterior samples. """
self._num_samples = thetas.shape[0]

wrapped = wrapper(self, 'model_sum')
for i in range(self._num_samples):
wrapped(None, thetas[i,:])

def __next__(self):
""" Update posterior expected model given the updated signal.
.. note::
You cannot make an iterator from an instance of this class.
"""
try:
self._model_sum
except AttributeError:
self._model_sum = self._signal.expected_counts.copy()
else:
self._model_sum = self._signal.expected_counts

@property
def model_sum(self):
""" Get the current posterior sum of the count numbers. """
return self._model_sum

@model_sum.deleter
def model_sum(self):
del self._model_sum

@property
def expected_counts(self):
""" Get the estimated posterior expectation of the count numbers. """
return self._model_sum / self._num_samples

@make_verbose('ResidualPlot object finalizing',
'ResidualPlot object finalized')
def finalize(self):
self._add_1d_residuals()

def _add_1d_residuals(self):
""" Display the residuals in the third panel. """

# Get the residuals
self._residuals = self.expected_counts - self._signal.data.counts
self._residuals /= _np.sqrt(self.expected_counts)
self._residuals = self._residuals.flatten()

# Plot them
self._ax_1d_residuals.hist(self._residuals, bins=self._nbins, density=True, alpha=0.5)

# Add the fit if needed
if self._plot_fit:

mean,std=norm.fit( self._residuals )
xmin, xmax = _np.min( self._residuals ) , _np.max( self._residuals )
x = _np.linspace(xmin, xmax, 100)
y = norm.pdf(x, mean, std)
self._ax_1d_residuals.plot( x , y , label=r'Fitted Normal $\sigma=$' f'{std:.3f}', color='darkblue' )
self._ax_1d_residuals.axvline(mean, color='k', linestyle='dashed', linewidth=2, label=r'$\mu=$' f'{mean:.5f}')
self._ax_1d_residuals.legend()


# Do the nice plotting
self._ax_1d_residuals.set_title('Residual histogram')
self._ax_1d_residuals.set_xlabel(r'$(c_{ik}-d_{ik})/\sqrt{c_{ik}}$')
self._ax_1d_residuals.set_ylabel('Density')
4 changes: 3 additions & 1 deletion xpsi/PostProcessing/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,11 60,13 @@
"SignalPlotter",
"PulsePlot",
"SpectrumPlot",
"ResidualPlot"]
"ResidualPlot",
"Residual1DPlot"]

from ._runs import Runs
from ._signalplotter import SignalPlotter
from ._residual import ResidualPlot
from ._1d_residual import Residual1DPlot
from ._pulse import PulsePlot
from ._spectrum import SpectrumPlot
from ._backends import NestedBackend
Expand Down
4 changes: 4 additions & 0 deletions xpsi/PostProcessing/_pulse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 118,7 @@ def __init__(self,
plot_truth=False,
truth_line_kwargs=None,
comp_truth_line_kwargs=None,
parameters_vector=None,
**kwargs):
super(PulsePlot, self).__init__(**kwargs)

Expand Down Expand Up @@ -201,6 202,9 @@ def __init__(self,
else:
self._expectation_line_kwargs = expectation_line_kwargs

if not parameters_vector is None:
self._parameters_vector = parameters_vector

if comp_expectation_line_kwargs is not None:
self._comp_expectation_line_kwargs = comp_expectation_line_kwargs
else:
Expand Down
Loading

0 comments on commit f3a7e81

Please sign in to comment.