Skip to content

Commit

Permalink
one dataset no list
Browse files Browse the repository at this point in the history
  • Loading branch information
@StFroese
StFroese committed Nov 28, 2024
1 parent e08888c commit cd8fc04
Show file tree
Hide file tree
Showing 6 changed files with 744 additions and 47 deletions.
Binary file added .DS_Store
View file
Binary file not shown.
616 changes: 616 additions & 0 deletions Test SciPy chi2.ipynb
View file

Large diffs are not rendered by default.

Binary file added examples/test.fits
View file
Binary file not shown.
54 changes: 54 additions & 0 deletions test.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
import astropy.units as u
from astropy import constants as const
from astropy.coordinates import SkyCoord
from gammapy.astro.darkmatter import JFactory, profiles
from gammapy.datasets import MapDataset
from gammapy.maps import WcsGeom, WcsNDMap

from titrate.upperlimits import ULFactory

if __name__ == "__main__":
filename = "$GAMMAPY_DATA/cta-1dc-gc/cta-1dc-gc.fits.gz"
dataset = MapDataset.read(filename, name="cta-dataset")

def geometry2d():
geom = WcsGeom.create(
skydir=SkyCoord(0, 0, unit="deg", frame="galactic"),
width=(4, 4),
binsz=0.04,
frame="galactic",
)

return geom

def ursa_major_ii_profile():
rhos = (
10**-1.1331
* const.M_sun.to(u.GeV, equivalencies=u.mass_energy())
/ u.pc**3
)
rs = 10**3.6317 * u.pc

profile = profiles.NFWProfile(r_s=rs, rho_s=rhos)
profile.DISTANCE_GC = 32 * u.kpc

return profile

jfactory = JFactory(
geom=geometry2d(),
profile=ursa_major_ii_profile(),
distance=ursa_major_ii_profile().DISTANCE_GC,
)
jfactor = jfactory.compute_differential_jfactor()
jfactor_map = WcsNDMap(geom=geometry2d(), data=jfactor.value, unit=jfactor.unit)

ul_factory = ULFactory(
dataset,
["b", "W", "tau", "mu"],
0.1 * u.TeV,
100 * u.TeV,
20,
jfactor_map,
max_workers=8,
)
ul_factory.compute()
64 changes: 31 additions & 33 deletions titrate/statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,6 @@
from gammapy.modeling import Fit
from scipy.stats import kstwo, norm

from titrate.datasets import AsimovMapDataset


class POIError(IndexError):
"""Parameter of interest is not defined in model"""
Expand Down Expand Up @@ -105,11 +103,11 @@ def check_for_pois(self):
return pois

def sigma(self):
if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to calculate"
" `sigma`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to calculate"
# " `sigma`"
# )
return np.sqrt(self.fit_result.covariance_result.matrix[0, 0])

def asympotic_approximation_pdf(
Expand All @@ -122,11 +120,11 @@ def asympotic_approximation_pdf(
* np.exp(-0.5 * (np.sqrt(ts_val)) ** 2)
)

if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to use the"
" `asympotic_approximation`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to use the"
# " `asympotic_approximation`"
# )
return (
1
/ (2 * np.sqrt(2 * np.pi * ts_val))
Expand All @@ -141,11 +139,11 @@ def asympotic_approximation_cdf(
if same:
return norm.cdf(np.sqrt(ts_val))

if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to use the"
" `asympotic_approximation`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to use the"
# " `asympotic_approximation`"
# )

return norm.cdf(np.sqrt(ts_val) - (poi_val - poi_true_val) / self.sigma())

Expand Down Expand Up @@ -186,7 +184,7 @@ def __init__(self, dataset, poi_name=""):
)

self.fit = Fit()
self.fit_result = self.fit.run(datasets=[self.dataset])
self.fit_result = self.fit.run(datasets=self.dataset)
self.poi_best = self.fit_result.parameters[self.poi_name].value
if self.poi_best < 0:
self.dataset.models.parameters[self.poi_name].scan_values = [0]
Expand Down Expand Up @@ -238,21 +236,21 @@ def check_for_pois(self):
return pois

def sigma(self):
if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to calculate"
" `sigma`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to calculate"
# " `sigma`"
# )
return np.sqrt(self.fit_result.covariance_result.matrix[0, 0])

def asympotic_approximation_pdf(
self, ts_val, poi_val, same=True, poi_true_val=None
):
if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to use the"
" `asympotic_approximation`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to use the"
# " `asympotic_approximation`"
# )

sigma = self.sigma()

Expand Down Expand Up @@ -287,11 +285,11 @@ def asympotic_approximation_pdf(
def asympotic_approximation_cdf(
self, ts_val, poi_val, same=True, poi_true_val=None
):
if not isinstance(self.dataset, AsimovMapDataset):
raise AsimovApproximationError(
"`dataset` must be an `AsimovMapDataset` in order to use the"
" `asympotic_approximation`"
)
# if not isinstance(self.dataset, AsimovMapDataset):
# raise AsimovApproximationError(
# "`dataset` must be an `AsimovMapDataset` in order to use the"
# " `asympotic_approximation`"
# )

sigma = self.sigma()
if same:
Expand Down
57 changes: 43 additions & 14 deletions titrate/upperlimits.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
import numpy as np
from astropy.table import QTable
from gammapy.astro.darkmatter import DarkMatterAnnihilationSpectralModel
from gammapy.datasets import Datasets
from gammapy.modeling import Fit
from gammapy.modeling.models import (
FoVBackgroundModel,
Expand Down Expand Up @@ -32,6 +31,7 @@ def __init__(
poi_name="scale",
cl=0.95,
cl_type="s",
analysis="3d",
):
self.measurement_dataset = measurement_dataset
self.poi_name = poi_name
Expand All @@ -47,6 +47,7 @@ def __init__(
raise ValueError('cl_type must be either "s+b" or "s"')
self.cl_type = cl_type
self.cl = cl
self.analysis = analysis

def compute(self):
poi_best = self.statistic.poi_best
Expand All @@ -55,35 +56,53 @@ def compute(self):
poi_ul = 1e-2
else:
poi_ul = poi_best + 0.01 * poi_best
while self.pvalue(poi_ul, cl_type=self.cl_type) > 1 - self.cl or np.isnan(
self.pvalue(poi_ul, cl_type=self.cl_type)
):
prev_pval = 0
while (
(pval := self.pvalue(poi_ul, cl_type=self.cl_type)) > 1 - self.cl
) or np.isnan(pval):
prev_pval = pval
poi_ul *= 2

if prev_pval == 0:
return np.nan
interp_ul_points = np.linspace(poi_ul / 2, poi_ul, 10)
interp_pvalues = np.array(
[
self.pvalue(interp_ul, cl_type=self.cl_type)
for interp_ul in interp_ul_points
]
).ravel()
interp_pvalues[0] = prev_pval
print(prev_pval, pval)
interpolation = interp1d(interp_ul_points, interp_pvalues - 1 + self.cl)
if interpolation(interp_ul_points[0]) < 0:
print(
interpolation(interp_ul_points[0]),
interpolation(interp_ul_points[-1]),
self.measurement_dataset,
)

poi_ul = brentq(interpolation, poi_ul / 2, poi_ul)

return poi_ul

def pvalue(self, poi_ul, cl_type):
pval_bkg = 0
if self.statistic.__class__.__name__ == "QTildeMuTestStatistic":
asimov_dataset = AsimovMapDataset.from_MapDataset(self.measurement_dataset)
if self.analysis == "3d":
asimov_dataset = AsimovMapDataset.from_MapDataset(
self.measurement_dataset
)
elif self.analysis == "1d":
asimov_dataset = AsimovSpectralDataset.from_SpectralDataset(
self.measurement_dataset
)
asimov_dataset.models.parameters[self.poi_name].value = poi_ul
asimov_dataset.fake()
if self.analysis == "3d":
asimov_dataset.fake()
if self.analysis == "1d":
asimov_dataset.fake()

statistic = self.statistic.__class__(asimov_dataset, poi_name=self.poi_name)
ts_val, valid = self.statistic.evaluate(
poi_ul
Expand All @@ -94,11 +113,21 @@ def pvalue(self, poi_ul, cl_type):
pval_sig_bkg = statistic.pvalue(poi_ul, ts_val=ts_val)

if cl_type == "s":
no_signal_asimov_dataset = AsimovMapDataset.from_MapDataset(
self.measurement_dataset
)
if self.analysis == "3d":
no_signal_asimov_dataset = AsimovMapDataset.from_MapDataset(
self.measurement_dataset
)
elif self.analysis == "1d":
no_signal_asimov_dataset = (
AsimovSpectralDataset.from_SpectralDataset(
self.measurement_dataset
)
)
no_signal_asimov_dataset.models.parameters[self.poi_name].value = 0
no_signal_asimov_dataset.fake()
if self.analysis == "3d":
no_signal_asimov_dataset.fake()
if self.analysis == "1d":
no_signal_asimov_dataset.fake()
no_signal_statistic = self.statistic.__class__(
no_signal_asimov_dataset, poi_name=self.poi_name
)
Expand All @@ -121,10 +150,9 @@ def pvalue(self, poi_ul, cl_type):
def expected_uls(self):
# Create asimov dataset
# asimov_dataset = AsimovMapDataset.from_MapDataset(self.measurement_dataset)
asimov_dataset = []
for dataset in self.measurement_dataset:
asimov_dataset.append(AsimovSpectralDataset.from_SpectralDataset(dataset))
asimov_dataset = Datasets(asimov_dataset)
asimov_dataset = AsimovSpectralDataset.from_SpectralDataset(
self.measurement_dataset
)
# asimov_dataset =
# AsimovSpectralDataset.from_SpectralDatasets(self.measurement_dataset)
asimov_dataset.models.parameters[self.poi_name].value = 0
Expand Down Expand Up @@ -183,6 +211,7 @@ def __init__(
self.max_workers = max_workers
self.kwargs["cl"] = self.cl
self.kwargs["cl_type"] = self.cl_type
self.kwargs["analysis"] = self.analysis
self.uls = None
self.expected_uls = None

Expand Down

0 comments on commit cd8fc04

Please sign in to comment.