Skip to content

Commit b76e82d

Browse files
marialiubarskamarialiubarska
authored
Merge pull request #766 from icecube/covar
Take into account covariance matrix for daemonflux parameters
2 parents 480c067 + b8ec5fb commit b76e82d

File tree

6 files changed

+546
-446
lines changed

6 files changed

+546
-446
lines changed

‎pisa/core/param.py

+19-2
Original file line numberDiff line numberDiff line change
@@ -1380,8 +1380,25 @@ def priors_penalty(self, metric):
13801380
penalty : float sum of all parameters' prior values
13811381
13821382
"""
1383-
return np.sum([obj.prior_penalty(metric=metric)
1384-
for obj in self._params])
1383+
1384+
# if daemonflux stage is not used use std priors penalty
1385+
if not "daemon_chi2" in self.names:
1386+
priors_sum = np.sum([obj.prior_penalty(metric=metric)
1387+
for obj in self._params])
1388+
1389+
# else switch daemon flux params penalty to the one drom daemonflux
1390+
# (which takes into account covariance)
1391+
else:
1392+
# normal (non-correlated) penalty for non-daemonflux params
1393+
priors_sum = np.sum([obj.prior_penalty(metric=metric)
1394+
for obj in self._params if "daemon_" not in obj.name])
1395+
1396+
# conversion factor between chi2 and llh
1397+
conv_factor = -0.5 if metric in LLH_METRICS else 1
1398+
# add daemonflux calcualted chi2 penalty
1399+
priors_sum += conv_factor * self._by_name["daemon_chi2"].value.m_as("dimensionless")
1400+
1401+
return priors_sum
13851402

13861403
def priors_penalties(self, metric):
13871404
"""Return the prior penalties for each param at their current values.

‎pisa/core/pipeline.py

+7
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,13 @@ def __init__(self, config, profile=False):
113113
self._init_stages()
114114
self._source_code_hash = None
115115

116+
# check in case someone decided to add a non-daemonflux parameter with daemon_
117+
# in it, which would potentially make penalty calculation incorrect
118+
if "daemon_chi2" in self.params.names:
119+
num_daemon_params = len([name for name in self.params.names if "daemon_" in name])
120+
assert num_daemon_params == self.params["daemon_params_len"].value.m_as("dimensionless"), \
121+
'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!'
122+
116123
self._covariance_set = False
117124

118125
def __repr__(self):

‎pisa/stages/flux/daemon_flux.py

+63-83
Original file line numberDiff line numberDiff line change
@@ -1,78 +1,81 @@
11
"""
2-
Implementation of DEAMON flux (https://arxiv.org/abs/2303.00022)
2+
Implementation of DAEMON flux (https://arxiv.org/abs/2303.00022)
33
by Juan Pablo Yañez and Anatoli Fedynitch for use in PISA.
44
55
Maria Liubarska, J.P. Yanez 2023
66
"""
77

88
import numpy as np
99
from daemonflux import Flux
10+
from daemonflux import __version__ as daemon_version
1011

1112
from pisa import FTYPE, TARGET
1213
from pisa.core.stage import Stage
14+
from pisa.core.param import Param
1315
from pisa.utils.log import logging
1416
from pisa.utils.profiler import profile
1517
from numba import jit
1618
from scipy import interpolate
19+
from packaging.version import Version
1720

1821

1922
class daemon_flux(Stage):
2023
"""
21-
DEAMON fulx stage
24+
DAEMON flux stage
2225
2326
Parameters
2427
----------
2528
2629
params: ParamSet
2730
Must have parameters: .. ::
2831
29-
K_158G : quantity (dimensionless)
32+
daemon_K_158G : quantity (dimensionless)
3033
31-
K_2P : quantity (dimensionless)
34+
daemon_K_2P : quantity (dimensionless)
3235
33-
K_31G : quantity (dimensionless)
36+
daemon_K_31G : quantity (dimensionless)
3437
35-
antiK_158G : quantity (dimensionless)
38+
daemon_antiK_158G : quantity (dimensionless)
3639
37-
antiK_2P : quantity (dimensionless)
40+
daemon_antiK_2P : quantity (dimensionless)
3841
39-
antiK_31G : quantity (dimensionless)
42+
daemon_antiK_31G : quantity (dimensionless)
4043
41-
n_158G : quantity (dimensionless)
44+
daemon_n_158G : quantity (dimensionless)
4245
43-
n_2P : quantity (dimensionless)
46+
daemon_n_2P : quantity (dimensionless)
4447
45-
p_158G : quantity (dimensionless)
48+
daemon_p_158G : quantity (dimensionless)
4649
47-
p_2P : quantity (dimensionless)
50+
daemon_p_2P : quantity (dimensionless)
4851
49-
pi_158G : quantity (dimensionless)
52+
daemon_pi_158G : quantity (dimensionless)
5053
51-
pi_20T : quantity (dimensionless)
54+
daemon_pi_20T : quantity (dimensionless)
5255
53-
pi_2P : quantity (dimensionless)
56+
daemon_pi_2P : quantity (dimensionless)
5457
55-
pi_31G : quantity (dimensionless)
58+
daemon_pi_31G : quantity (dimensionless)
5659
57-
antipi_158G : quantity (dimensionless)
60+
daemon_antipi_158G : quantity (dimensionless)
5861
59-
antipi_20T : quantity (dimensionless)
62+
daemon_antipi_20T : quantity (dimensionless)
6063
61-
antipi_2P : quantity (dimensionless)
64+
daemon_antipi_2P : quantity (dimensionless)
6265
63-
antipi_31G : quantity (dimensionless)
66+
daemon_antipi_31G : quantity (dimensionless)
6467
65-
GSF_1 : quantity (dimensionless)
68+
daemon_GSF_1 : quantity (dimensionless)
6669
67-
GSF_2 : quantity (dimensionless)
70+
daemon_GSF_2 : quantity (dimensionless)
6871
69-
GSF_3 : quantity (dimensionless)
72+
daemon_GSF_3 : quantity (dimensionless)
7073
71-
GSF_4 : quantity (dimensionless)
74+
daemon_GSF_4 : quantity (dimensionless)
7275
73-
GSF_5 : quantity (dimensionless)
76+
daemon_GSF_5 : quantity (dimensionless)
7477
75-
GSF_6 : quantity (dimensionless)
78+
daemon_GSF_6 : quantity (dimensionless)
7679
7780
"""
7881

@@ -81,70 +84,43 @@ def __init__(
8184
**std_kwargs,
8285
):
8386

84-
self.deamon_params = ['K_158G',
85-
'K_2P',
86-
'K_31G',
87-
'antiK_158G',
88-
'antiK_2P',
89-
'antiK_31G',
90-
'n_158G',
91-
'n_2P',
92-
'p_158G',
93-
'p_2P',
94-
'pi_158G',
95-
'pi_20T',
96-
'pi_2P',
97-
'pi_31G',
98-
'antipi_158G',
99-
'antipi_20T',
100-
'antipi_2P',
101-
'antipi_31G',
102-
'GSF_1',
103-
'GSF_2',
104-
'GSF_3',
105-
'GSF_4',
106-
'GSF_5',
107-
'GSF_6',
108-
]
109-
110-
self.deamon_names = ['K+_158G',
111-
'K+_2P',
112-
'K+_31G',
113-
'K-_158G',
114-
'K-_2P',
115-
'K-_31G',
116-
'n_158G',
117-
'n_2P',
118-
'p_158G',
119-
'p_2P',
120-
'pi+_158G',
121-
'pi+_20T',
122-
'pi+_2P',
123-
'pi+_31G',
124-
'pi-_158G',
125-
'pi-_20T',
126-
'pi-_2P',
127-
'pi-_31G',
128-
'GSF_1',
129-
'GSF_2',
130-
'GSF_3',
131-
'GSF_4',
132-
'GSF_5',
133-
'GSF_6',
134-
]
87+
# first have to check daemonflux package version is >=0.8.0
88+
# (have to do this to ensure chi2 prior penalty is calculated correctly)
89+
if Version(daemon_version) < Version("0.8.0"):
90+
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'")
91+
raise Exception('detected daemonflux version < 0.8.0')
92+
93+
# create daemonflux Flux object
94+
self.flux_obj = Flux(location='IceCube', use_calibration=True)
95+
96+
# get parameter names from daemonflux
97+
self.daemon_names = self.flux_obj.params.known_parameters
98+
99+
# make parameter names pisa config compatible and add prefix
100+
self.daemon_params = ['daemon_'+p.replace('pi+','pi').replace('pi-','antipi').replace('K+','K').replace('K-','antiK')
101+
for p in self.daemon_names]
102+
103+
# add daemon_chi2 internal parameter to carry on chi2 penalty from daemonflux (using covar. matrix)
104+
daemon_chi2 = Param(name='daemon_chi2', nominal_value=0.,
105+
value=0., prior=None, range=None, is_fixed=True)
106+
107+
# saving number of parameters into a internal param in order to check that we don't have
108+
# non-daemonflux params with 'daemon_' in their name, which will make prior penalty calculation incorrect
109+
daemon_params_len = Param(name='daemon_params_len', nominal_value=len(self.daemon_names)+2,
110+
value=len(self.daemon_names)+2, prior=None, range=None, is_fixed=True)
111+
112+
std_kwargs['params'].update([daemon_chi2,daemon_params_len])
135113

136114
# init base class
137115
super(daemon_flux, self).__init__(
138-
expected_params=tuple(self.deamon_params),
116+
expected_params=tuple(self.daemon_params+['daemon_chi2','daemon_params_len']),
139117
**std_kwargs,
140118
)
141119

142120
def setup_function(self):
143121

144122
self.data.representation = self.calc_mode
145123

146-
self.flux_obj = Flux(location='IceCube')
147-
148124
for container in self.data:
149125
container['nu_flux'] = np.empty((container.size, 2), dtype=FTYPE)
150126

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

156132
# get modified parameters (in units of sigma)
157133
modif_param_dict = {}
158-
for i,k in enumerate(self.deamon_params):
159-
modif_param_dict[self.deamon_names[i]] = getattr(self.params, k).value.m_as("dimensionless")
134+
for i,k in enumerate(self.daemon_params):
135+
modif_param_dict[self.daemon_names[i]] = getattr(self.params, k).value.m_as("dimensionless")
136+
137+
# update chi2 parameter
138+
self.params['daemon_chi2'].value = self.flux_obj.chi2(modif_param_dict)
160139

140+
# compute flux maps
161141
flux_map_numu = make_2d_flux_map(self.flux_obj,
162142
particle = 'numu',
163143
params = modif_param_dict)
@@ -219,4 +199,4 @@ def evaluate_flux_map(flux_map, true_energy, true_coszen):
219199

220200
uconv = true_energy**-3 * 1e4
221201
return flux_map.ev(true_energy, true_coszen) * uconv
222-
202+

0 commit comments

Comments
 (0)