Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Take into account covariance matrix for daemonflux parameters #766

Merged
merged 9 commits into from
Apr 23, 2024
Merged

Take into account covariance matrix for daemonflux parameters #766

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a612062
add support for daemonflux penalty term calculation taking into accou…
Apr 11, 2024
ef2c0e5
fix list of daemonflux param names to follow param naming, not daemon…
marialiubarska Apr 11, 2024
c3d7750
update daemonflux example script
marialiubarska Apr 11, 2024
3803a82
add version condition for daemonflux and add packaging to required de…
marialiubarska Apr 11, 2024
defc16e
break if daemonflux version is not updated, in case people didn't rei…
marialiubarska Apr 11, 2024
e95ae41
fix misspelled names (cosmetic)
marialiubarska Apr 17, 2024
098d603
change the way we identify daemonflux params and also copy list of pa…
marialiubarska Apr 17, 2024
9abeaf0
break if someone creates a parameter with 'daemon_' in it's name whic…
marialiubarska Apr 23, 2024
b8ec5fb
add conversion of covariance prior penalty to llh for daemonflux and …
marialiubarska Apr 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 19 additions & 2 deletions pisa/core/param.py
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using this branch, I tried to calculate the pull penalty, and it seems to give a number as expected. However, I have the following observations/concerns:

  • There are two other functions also params.priors_chi2 and params.priors_llh. I believe their output are supposed be in agreement with params.priors_penalty(metric='mod_chi2'). However, those two still give the priors without incorporating the covariance matrix of daemonflux.

  • params.priors_penalty(metric='llh') gives output same as params.priors_penalty(metric='mod_chi2') for modification to daemonflux parameters. If I modify other parameters, let's say dom_eff then params.priors_penalty(metric='llh') give me a value which is half of params.priors_penalty(metric='mod_chi2') as expected. I think this can cause an issue if someone performs minimization using llh as the metric. Therefore, we need introduce a factor of 1/2 if the metric is llh.

  • If I modify the values of the daemonflux parameters in params then params.priors_chi2 and params.priors_llh give me updated priors immediately. But params.priors_penalty(metric='llh') gives me an updated value only after I run template_maker.get_outputs(return_sum=True). Otherwise old value is given. If someone loads a fitted paramset and tries to calculate the params.priors_penalty without generating a template using get_outputs then the contribution from daemonflux parameter via covariance matrix will be missing but other parameters will still give a contribution to the prior penalty. This can cause confusion and go unnoticed leading to a wrong prior penalty.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JanWeldert Calculation is ok when the metric is chi2 but these inconsistencies mentioned in the previous message can cause confusion and issues. More specifically, if some try to use LLH, then they can get the wrong numbers.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, these are good points.

  • You could modify params.priors_chi2 and params.priors_llh in the same way as params.priors_penalty, but we could also think about removing the functions completely since they are just calling the prior penalty function of the Param class anyway. params.priors_penalties will also ignore the correlation, would be good to exclude the daemon params here too.
  • Here you can simply check if metric in LLH_METRICS: and multiply the chi2 value by 0.5 if so. You can also make sure that metric in CHI2_METRICS otherwise.
  • Not sure if I understand that point. If you want to update the daemonflux chi2. you have to run the distribution maker since it is calculated in the flux stage. The individual contributions of the daemonflux parameters are ignored by params.priors_penalty (as they should be) but currently taken into account by params.priors_chi2 and params.priors_llh because they don't know about the correlation (your first point).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@anilak41 Thanks for pointing out the issue with LLH conversion, I added a fix for this

Original file line number Diff line number Diff line change
Expand Up @@ -1380,8 +1380,25 @@ def priors_penalty(self, metric):
penalty : float sum of all parameters' prior values

"""
return np.sum([obj.prior_penalty(metric=metric)
for obj in self._params])

# if daemonflux stage is not used use std priors penalty
if not "daemon_chi2" in self.names:
priors_sum = np.sum([obj.prior_penalty(metric=metric)
for obj in self._params])

# else switch daemon flux params penalty to the one drom daemonflux
# (which takes into account covariance)
else:
# normal (non-correlated) penalty for non-daemonflux params
priors_sum = np.sum([obj.prior_penalty(metric=metric)
for obj in self._params if "daemon_" not in obj.name])

# conversion factor between chi2 and llh
conv_factor = -0.5 if metric in LLH_METRICS else 1
# add daemonflux calcualted chi2 penalty
priors_sum += conv_factor * self._by_name["daemon_chi2"].value.m_as("dimensionless")

return priors_sum

def priors_penalties(self, metric):
"""Return the prior penalties for each param at their current values.
Expand Down
7 changes: 7 additions & 0 deletions pisa/core/pipeline.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,13 @@ def __init__(self, config, profile=False):
self._init_stages()
self._source_code_hash = None

# check in case someone decided to add a non-daemonflux parameter with daemon_
# in it, which would potentially make penalty calculation incorrect
if "daemon_chi2" in self.params.names:
num_daemon_params = len([name for name in self.params.names if "daemon_" in name])
assert num_daemon_params == self.params["daemon_params_len"].value.m_as("dimensionless"), \
'Incorrect number of parameters with "daemon_" in their name detected. Non-daemonflux parameters can not have "daemon_" in their name. Rename your non-daemonflux parameters which do not comly!'

self._covariance_set = False

def __repr__(self):
Expand Down
146 changes: 63 additions & 83 deletions pisa/stages/flux/daemon_flux.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,78 +1,81 @@
"""
Implementation of DEAMON flux (https://arxiv.org/abs/2303.00022)
Implementation of DAEMON flux (https://arxiv.org/abs/2303.00022)
by Juan Pablo Yañez and Anatoli Fedynitch for use in PISA.

Maria Liubarska, J.P. Yanez 2023
"""

import numpy as np
from daemonflux import Flux
from daemonflux import __version__ as daemon_version

from pisa import FTYPE, TARGET
from pisa.core.stage import Stage
from pisa.core.param import Param
from pisa.utils.log import logging
from pisa.utils.profiler import profile
from numba import jit
from scipy import interpolate
from packaging.version import Version


class daemon_flux(Stage):
"""
DEAMON fulx stage
DAEMON flux stage

Parameters
----------

params: ParamSet
Must have parameters: .. ::

K_158G : quantity (dimensionless)
daemon_K_158G : quantity (dimensionless)

K_2P : quantity (dimensionless)
daemon_K_2P : quantity (dimensionless)

K_31G : quantity (dimensionless)
daemon_K_31G : quantity (dimensionless)

antiK_158G : quantity (dimensionless)
daemon_antiK_158G : quantity (dimensionless)

antiK_2P : quantity (dimensionless)
daemon_antiK_2P : quantity (dimensionless)

antiK_31G : quantity (dimensionless)
daemon_antiK_31G : quantity (dimensionless)

n_158G : quantity (dimensionless)
daemon_n_158G : quantity (dimensionless)

n_2P : quantity (dimensionless)
daemon_n_2P : quantity (dimensionless)

p_158G : quantity (dimensionless)
daemon_p_158G : quantity (dimensionless)

p_2P : quantity (dimensionless)
daemon_p_2P : quantity (dimensionless)

pi_158G : quantity (dimensionless)
daemon_pi_158G : quantity (dimensionless)

pi_20T : quantity (dimensionless)
daemon_pi_20T : quantity (dimensionless)

pi_2P : quantity (dimensionless)
daemon_pi_2P : quantity (dimensionless)

pi_31G : quantity (dimensionless)
daemon_pi_31G : quantity (dimensionless)

antipi_158G : quantity (dimensionless)
daemon_antipi_158G : quantity (dimensionless)

antipi_20T : quantity (dimensionless)
daemon_antipi_20T : quantity (dimensionless)

antipi_2P : quantity (dimensionless)
daemon_antipi_2P : quantity (dimensionless)

antipi_31G : quantity (dimensionless)
daemon_antipi_31G : quantity (dimensionless)

GSF_1 : quantity (dimensionless)
daemon_GSF_1 : quantity (dimensionless)

GSF_2 : quantity (dimensionless)
daemon_GSF_2 : quantity (dimensionless)

GSF_3 : quantity (dimensionless)
daemon_GSF_3 : quantity (dimensionless)

GSF_4 : quantity (dimensionless)
daemon_GSF_4 : quantity (dimensionless)

GSF_5 : quantity (dimensionless)
daemon_GSF_5 : quantity (dimensionless)

GSF_6 : quantity (dimensionless)
daemon_GSF_6 : quantity (dimensionless)

"""

Expand All @@ -81,70 +84,43 @@ def __init__(
**std_kwargs,
):

self.deamon_params = ['K_158G',
'K_2P',
'K_31G',
'antiK_158G',
'antiK_2P',
'antiK_31G',
'n_158G',
'n_2P',
'p_158G',
'p_2P',
'pi_158G',
'pi_20T',
'pi_2P',
'pi_31G',
'antipi_158G',
'antipi_20T',
'antipi_2P',
'antipi_31G',
'GSF_1',
'GSF_2',
'GSF_3',
'GSF_4',
'GSF_5',
'GSF_6',
]

self.deamon_names = ['K+_158G',
'K+_2P',
'K+_31G',
'K-_158G',
'K-_2P',
'K-_31G',
'n_158G',
'n_2P',
'p_158G',
'p_2P',
'pi+_158G',
'pi+_20T',
'pi+_2P',
'pi+_31G',
'pi-_158G',
'pi-_20T',
'pi-_2P',
'pi-_31G',
'GSF_1',
'GSF_2',
'GSF_3',
'GSF_4',
'GSF_5',
'GSF_6',
]
# first have to check daemonflux package version is >=0.8.0
# (have to do this to ensure chi2 prior penalty is calculated correctly)
if Version(daemon_version) < Version("0.8.0"):
logging.fatal("Detected daemonflux version below 0.8.0! This will lead to incorrect penalty calculation. You must update your daemonflux package to use this stage. You can do it by running 'pip install daemonflux --upgrade'")
raise Exception('detected daemonflux version < 0.8.0')

# create daemonflux Flux object
self.flux_obj = Flux(location='IceCube', use_calibration=True)

# get parameter names from daemonflux
self.daemon_names = self.flux_obj.params.known_parameters

# make parameter names pisa config compatible and add prefix
self.daemon_params = ['daemon_'+p.replace('pi+','pi').replace('pi-','antipi').replace('K+','K').replace('K-','antiK')
for p in self.daemon_names]

# add daemon_chi2 internal parameter to carry on chi2 penalty from daemonflux (using covar. matrix)
daemon_chi2 = Param(name='daemon_chi2', nominal_value=0.,
value=0., prior=None, range=None, is_fixed=True)

# saving number of parameters into a internal param in order to check that we don't have
# non-daemonflux params with 'daemon_' in their name, which will make prior penalty calculation incorrect
daemon_params_len = Param(name='daemon_params_len', nominal_value=len(self.daemon_names)+2,
value=len(self.daemon_names)+2, prior=None, range=None, is_fixed=True)

std_kwargs['params'].update([daemon_chi2,daemon_params_len])

# init base class
super(daemon_flux, self).__init__(
expected_params=tuple(self.deamon_params),
expected_params=tuple(self.daemon_params+['daemon_chi2','daemon_params_len']),
**std_kwargs,
)

def setup_function(self):

self.data.representation = self.calc_mode

self.flux_obj = Flux(location='IceCube')

for container in self.data:
container['nu_flux'] = np.empty((container.size, 2), dtype=FTYPE)

Expand All @@ -155,9 +131,13 @@ def compute_function(self):

# get modified parameters (in units of sigma)
modif_param_dict = {}
for i,k in enumerate(self.deamon_params):
modif_param_dict[self.deamon_names[i]] = getattr(self.params, k).value.m_as("dimensionless")
for i,k in enumerate(self.daemon_params):
modif_param_dict[self.daemon_names[i]] = getattr(self.params, k).value.m_as("dimensionless")

# update chi2 parameter
self.params['daemon_chi2'].value = self.flux_obj.chi2(modif_param_dict)

# compute flux maps
flux_map_numu = make_2d_flux_map(self.flux_obj,
particle = 'numu',
params = modif_param_dict)
Expand Down Expand Up @@ -219,4 +199,4 @@ def evaluate_flux_map(flux_map, true_energy, true_coszen):

uconv = true_energy**-3 * 1e4
return flux_map.ev(true_energy, true_coszen) * uconv


Loading
Loading