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 5 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
45 changes: 43 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
Collaborator

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 @@ -862,6 +862,33 @@ def __init__(self, *args):
# I think because the changed units are cached in the object (Philipp)
self.normalize_values = True

# store list of daemonflux params names for prior penalty calculation
self._daemon_names = ['K_158G',
Copy link

Choose a reason for hiding this comment

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

To ensure future compatibility it would be good read the daemonflux parameter names from Flux().params.known_parameters instead of a hard coded list. If, for some reason, you want to keep a fixed parameter list, you may want to implement a check if what is returned by Flux().params.known_parameters is the same as what you expect it to be.

Copy link

Choose a reason for hiding this comment

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

to replace the names use simple string ops:

self._daemon_names = [
     p.replace("pi+","pi").replace("pi-","antipi").replace(
          "K+","K").replace("K-","antiK")
      for p in Flux().params.known_parameters
]```

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The idea was to limit daemonflux dependency to the corresponding stage, that's why I wanted to avoid calling it directly in ParamSet class (that's also why chi2 value is stored into a parameter inside the stage). But I also had to add daemonflux to required packages anyway for automatic installation to work and now to restrict the version, so maybe there is no harm in doing it this way.
I will look at this tomorrow.

I didn't think about renaming the parameters automatically, that's a good suggestion, thank you!

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Maybe another solution which avoids storing the parameter list altogether is to just add a prefix to all daemonflux parameters names, since we have to rename them anyway to remove symbols, then I can just look for that prefix as a condition to skip these parameters...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

maybe I should just do that

Copy link
Collaborator

Choose a reason for hiding this comment

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

Yeah, not having a hard coded list in param.py would be better

'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',
]

@property
def has_derived(self)->bool:
"""
Expand Down Expand Up @@ -1380,8 +1407,22 @@ 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 obj.name not in self._daemon_names])
# add daemonflux calcualted chi2 penalty
priors_sum += 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
24 changes: 21 additions & 3 deletions pisa/stages/flux/daemon_flux.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,16 @@

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):
Expand Down Expand Up @@ -81,6 +84,12 @@ def __init__(
**std_kwargs,
):

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

self.deamon_params = ['K_158G',
'K_2P',
'K_31G',
Expand Down Expand Up @@ -133,17 +142,22 @@ def __init__(
'GSF_6',
]

# 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)
std_kwargs['params'].update(daemon_chi2)

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

def setup_function(self):

self.data.representation = self.calc_mode

self.flux_obj = Flux(location='IceCube')
self.flux_obj = Flux(location='IceCube', use_calibration=True)

for container in self.data:
container['nu_flux'] = np.empty((container.size, 2), dtype=FTYPE)
Expand All @@ -158,6 +172,10 @@ def compute_function(self):
for i,k in enumerate(self.deamon_params):
modif_param_dict[self.deamon_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 +237,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


474 changes: 261 additions & 213 deletions pisa_examples/test_daemonflux_stage.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

3 changes: 2 additions & 1 deletion setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,8 @@
'scikit-learn<=1.1.2',
'pyarrow',
'tqdm',
'daemonflux',
'daemonflux>=0.8.0',
'packaging',
]

EXTRAS_REQUIRE = {
Expand Down
Loading