gammapy__gammapy-5674

{
  "status": "ok",
  "stdout": "gammapy/estimators/points/:\n__init__.py\ncore.py\nlightcurve.py\nprofile.py\nsed.py\nsensitivity.py\ntests\n\ngammapy/estimators/points/tests:\n__init__.py\ntest_core.py\ntest_lightcurve.py\ntest_profile.py\ntest_sed.py\ntest_sensitivity.py\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport logging\nfrom copy import deepcopy\nimport numpy as np\nfrom scipy import stats\nfrom scipy.interpolate import interp1d\nfrom scipy.optimize import minimize\nfrom astropy.io import fits\nfrom astropy.io.registry import IORegistryError\nfrom astropy.table import Table, vstack\nfrom astropy.time import Time\nfrom astropy.visualization import quantity_support\nimport matplotlib.pyplot as plt\nfrom gammapy.data import GTI\nfrom gammapy.maps import Map, MapAxis, Maps, RegionNDMap, TimeMapAxis\nfrom gammapy.maps.axes import UNIT_STRING_FORMAT, flat_if_equal\nfrom gammapy.modeling.models import TemplateSpectralModel\nfrom gammapy.modeling.models.spectral import scale_plot_flux\nfrom gammapy.modeling.scipy import stat_profile_ul_scipy\nfrom gammapy.utils.interpolation import interpolate_profile\nfrom gammapy.utils.scripts import make_path\nfrom gammapy.utils.table import table_standardise_units_copy\nfrom gammapy.utils.time import time_ref_to_dict\nfrom ..map.core import DEFAULT_UNIT, FluxMaps\n\n__all__ = [\"FluxPoints\"]\n\nlog = logging.getLogger(__name__)\n\n\ndef squash_fluxpoints(flux_point, axis):\n    \"\"\"Squash a `~FluxPoints` object into one point.\n\n    The log-likelihoods profiles in each bin are summed\n    to compute the resultant quantities. Stat profiles\n    must be present on the `~FluxPoints` object for\n    this method to work.\n    \"\"\"\n    value_scan = flux_point.stat_scan.geom.axes[\"norm\"].center\n    stat_scan = np.sum(flux_point.stat_scan.data, axis=0).ravel()\n    f = interp1d(value_scan, stat_scan, kind=\"quadratic\", bounds_error=False)\n    f = interpolate_profile(value_scan, stat_scan)\n    minimizer = minimize(\n        f,\n        x0=value_scan[int(len(value_scan) / 2)],\n        bounds=[(value_scan[0], value_scan[-1])],\n        method=\"L-BFGS-B\",\n    )\n\n    maps = Maps()\n    geom = flux_point.geom.to_image()\n    if axis.name != \"energy\":\n        geom = geom.to_cube([flux_point.geom.axes[\"energy\"]])\n\n    maps[\"norm\"] = Map.from_geom(geom, data=minimizer.x)\n    maps[\"norm_err\"] = Map.from_geom(\n        geom, data=np.sqrt(minimizer.hess_inv.todense()).reshape(geom.data_shape)\n    )\n    maps[\"n_dof\"] = Map.from_geom(geom, data=flux_point.geom.axes[axis.name].nbin)\n\n    if \"norm_ul\" in flux_point.available_quantities:\n        delta_ts = flux_point.meta.get(\"n_sigma_ul\", 2) ** 2\n        ul = stat_profile_ul_scipy(value_scan, stat_scan, delta_ts=delta_ts)\n        maps[\"norm_ul\"] = Map.from_geom(geom, data=ul.value)\n\n    maps[\"stat\"] = Map.from_geom(geom, data=f(minimizer.x))\n\n    geom_scan = geom.to_cube([MapAxis.from_nodes(value_scan, name=\"norm\")])\n    maps[\"stat_scan\"] = Map.from_geom(\n        geom=geom_scan, data=stat_scan.reshape(geom_scan.data_shape)\n    )\n    try:\n        maps[\"stat_null\"] = Map.from_geom(geom, data=np.sum(flux_point.stat_null.data))\n        maps[\"ts\"] = maps[\"stat_null\"] - maps[\"stat\"]\n    except AttributeError:\n        log.info(\n            \"Stat null info not present on original FluxPoints object. TS not computed\"\n        )\n\n    maps[\"success\"] = Map.from_geom(geom=geom, data=minimizer.success, dtype=bool)\n\n    combined_fp = FluxPoints.from_maps(\n        maps=maps,\n        sed_type=flux_point.sed_type_init,\n        reference_model=flux_point.reference_model,\n        gti=flux_point.gti,\n        meta=flux_point.meta,\n    )\n    return combined_fp\n\n\nclass FluxPoints(FluxMaps):\n    \"\"\"Flux points container.\n\n    The supported formats are described here: :ref:`gadf:flux-points`.\n\n    In summary, the following formats and minimum required columns are:\n\n    * Format ``dnde``: columns ``e_ref`` and ``dnde``\n    * Format ``e2dnde``: columns ``e_ref``, ``e2dnde``",
  "start_line": 1,
  "end_line": 100,
  "total_lines": 872,
  "is_truncated": false
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "        Parameters\n        ----------\n        axis_new : `MapAxis` or `TimeMapAxis`\n            The new axis to resample along\n\n        Returns\n        -------\n        flux_points : `~gammapy.estimators.FluxPoints`\n            A new FluxPoints object with modified axis.\n        \"\"\"\n        if not self.has_stat_profiles:\n            raise ValueError(\"Stat profiles not present, rebinning is not possible\")\n\n        fluxpoints = []\n        for edge_min, edge_max in zip(axis_new.edges_min, axis_new.edges_max):\n            if isinstance(axis_new, TimeMapAxis):\n                edge_min = edge_min + axis_new.reference_time\n                edge_max = edge_max + axis_new.reference_time\n            fp = self.slice_by_coord({axis_new.name: slice(edge_min, edge_max)})\n            fp_new = squash_fluxpoints(fp, axis_new)\n            fluxpoints.append(fp_new)\n\n        return self.__class__.from_stack(fluxpoints, axis=axis_new)",
  "start_line": 850,
  "end_line": 872,
  "total_lines": 872,
  "is_truncated": false
}

{
  "status": "ok",
  "stdout": "gammapy/estimators/points/tests/test_core.py\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/tests/test_core.py",
  "content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport pytest\nimport numpy as np\nfrom numpy.testing import assert_allclose\nimport astropy.units as u\nfrom astropy.table import Table\nfrom astropy.time import Time\nimport matplotlib.pyplot as plt\nfrom gammapy.catalog.fermi import SourceCatalog3FGL, SourceCatalog4FGL\nfrom gammapy.data import GTI\nfrom gammapy.estimators import FluxPoints\nfrom gammapy.estimators.map.core import DEFAULT_UNIT\nfrom gammapy.estimators.utils import get_rebinned_axis\nfrom gammapy.maps import MapAxis\nfrom gammapy.modeling.models import PowerLawSpectralModel, SpectralModel\nfrom gammapy.utils.scripts import make_path\nfrom gammapy.utils.testing import (\n    assert_quantity_allclose,\n    mpl_plot_check,\n    requires_data,\n)\n\nFLUX_POINTS_FILES = [\n    \"diff_flux_points.ecsv\",\n    \"diff_flux_points.fits\",\n    \"flux_points.ecsv\",\n    \"flux_points.fits\",\n]\n\n\nclass LWTestModel(SpectralModel):\n    @staticmethod\n    def evaluate(x):\n        return 1e4 * np.exp(-6 * x)\n\n    def integral(self, xmin, xmax, **kwargs):\n        return -1.0 / 6 * 1e4 * (np.exp(-6 * xmax) - np.exp(-6 * xmin))\n\n    def inverse(self, y):\n        return -1.0 / 6 * np.log(y * 1e-4)\n\n\nclass XSqrTestModel(SpectralModel):\n    @staticmethod\n    def evaluate(x):\n        return x**2\n\n    def integral(self, xmin, xmax, **kwargs):\n        return 1.0 / 3 * (xmax**3 - xmin**2)\n\n    def inverse(self, y):\n        return np.sqrt(y)\n\n\nclass ExpTestModel(SpectralModel):\n    @staticmethod\n    def evaluate(x):\n        return np.exp(x * u.Unit(\"1 / TeV\"))\n\n    def integral(self, xmin, xmax, **kwargs):\n        return np.exp(xmax * u.Unit(\"1 / TeV\")) - np.exp(xmin * u.Unit(\"1 / TeV\"))\n\n    def inverse(self, y):\n        return np.log(y * u.TeV) * u.TeV\n\n\ndef test_energy_ref_lafferty():\n    \"\"\"\n    Tests Lafferty & Wyatt x-point method.\n\n    Using input function g(x) = 10^4 exp(-6x) against\n    check values from paper Lafferty & Wyatt. Nucl. Instr. and Meth. in Phys.\n    Res. A 355 (1995) 541-547, p. 542 Table 1\n    \"\"\"\n    # These are the results from the paper\n    desired = np.array([0.048, 0.190, 0.428, 0.762])\n\n    model = LWTestModel()\n    energy_min = np.array([0.0, 0.1, 0.3, 0.6])\n    energy_max = np.array([0.1, 0.3, 0.6, 1.0])\n    actual = FluxPoints._energy_ref_lafferty(model, energy_min, energy_max)\n    assert_allclose(actual, desired, atol=1e-3)\n\n\n@pytest.mark.xfail\ndef test_dnde_from_flux():\n    \"\"\"Tests y-value normalization adjustment method.\"\"\"\n    table = Table()\n    table[\"e_min\"] = np.array([10, 20, 30, 40])\n    table[\"e_max\"] = np.array([20, 30, 40, 50])\n    table[\"flux\"] = np.array([42, 52, 62, 72])  # 'True' integral flux in this test bin\n\n    # Get values\n    model = XSqrTestModel()\n    table[\"e_ref\"] = FluxPoints._energy_ref_lafferty(\n        model, table[\"e_min\"], table[\"e_max\"]\n    )\n    dnde = FluxPoints.from_table(table, reference_model=model)\n\n    # Set up test case comparison\n    dnde_model = model(table[\"e_ref\"])\n\n    # Test comparison result\n    desired = model.integral(table[\"e_min\"], table[\"e_max\"])\n    # Test output result\n    actual = table[\"flux\"] * (dnde_model / dnde)\n    # Compare\n    assert_allclose(actual, desired, rtol=1e-6)\n\n\n@pytest.mark.xfail\n@pytest.mark.parametrize(\"method\", [\"table\", \"lafferty\", \"log_center\"])\ndef test_compute_flux_points_dnde_exp(method):\n    \"\"\"\n    Tests against analytical result or result from a powerlaw.\n    \"\"\"\n    model = ExpTestModel()\n\n    energy_min = [1.0, 10.0] * u.TeV\n    energy_max = [10.0, 100.0] * u.TeV\n\n    table = Table()\n    table.meta[\"SED_TYPE\"] = \"flux\"\n    table[\"e_min\"] = energy_min\n    table[\"e_max\"] = energy_max\n\n    flux = model.integral(energy_min, energy_max)\n    table[\"flux\"] = flux\n\n    if method == \"log_center\":\n        energy_ref = np.sqrt(energy_min * energy_max)\n    elif method == \"table\":\n        energy_ref = [2.0, 20.0] * u.TeV\n    elif method == \"lafferty\":\n        energy_ref = FluxPoints._energy_ref_lafferty(model, energy_min, energy_max)\n\n    table[\"e_ref\"] = energy_ref\n\n    result = FluxPoints.from_table(table, reference_model=model)\n\n    # Test energy\n    actual = result.energy_ref\n    assert_quantity_allclose(actual, energy_ref, rtol=1e-8)\n\n    # Test flux\n    actual = result.dnde\n    desired = model(energy_ref)\n    assert_quantity_allclose(actual, desired, rtol=1e-8)\n\n",
  "start_line": 1,
  "end_line": 150,
  "total_lines": 447,
  "is_truncated": true
}

{
  "status": "ok",
  "stdout": "498:            if self.has_stat_profiles:\n819:        if not self.has_stat_profiles:\n860:        if not self.has_stat_profiles:\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "\n        >>> from gammapy.estimators import FluxPoints\n        >>> fp = FluxPoints.read(\"$GAMMAPY_DATA/hawc_crab/HAWC19_flux_points.fits\")\n        >>> table = fp.to_table(sed_type=\"flux\", formatted=True)\n        >>> print(table[:2])\n        e_ref e_min e_max     flux      flux_err    flux_ul      ts    sqrt_ts is_ul\n         TeV   TeV   TeV  1 / (s cm2) 1 / (s cm2) 1 / (s cm2)\n        ----- ----- ----- ----------- ----------- ----------- -------- ------- -----\n        1.334 1.000 1.780   1.423e-11   3.135e-13         nan 2734.000  52.288 False\n        2.372 1.780 3.160   5.780e-12   1.082e-13         nan 4112.000  64.125 False\n        \"\"\"\n        if sed_type is None:\n            sed_type = self.sed_type_init\n\n        if format is None:\n            format = self._guess_format()\n            log.info(\"Inferred format: \" + format)\n\n        if format == \"gadf-sed\":\n            # TODO: what to do with GTI info?\n            if not self.geom.axes.names == [\"energy\"]:\n                raise ValueError(\n                    \"Only flux points with a single energy axis \"\n                    \"can be converted to 'gadf-sed'\"\n                )\n\n            idx = (Ellipsis, 0, 0)\n            table = self.energy_axis.to_table(format=\"gadf-sed\")\n            table.meta[\"SED_TYPE\"] = sed_type\n\n            if not self.is_convertible_to_flux_sed_type:\n                table.remove_columns([\"e_min\", \"e_max\"])\n\n            if self.n_sigma_ul:\n                table.meta[\"UL_CONF\"] = np.round(\n                    1 - 2 * stats.norm.sf(self.n_sigma_ul), 7\n                )\n\n            if sed_type == \"likelihood\":\n                table[\"ref_dnde\"] = self.dnde_ref[idx]\n                table[\"ref_flux\"] = self.flux_ref[idx]\n                table[\"ref_eflux\"] = self.eflux_ref[idx]\n\n            for quantity in self.all_quantities(sed_type=sed_type):\n                data = getattr(self, quantity, None)\n                if data:\n                    table[quantity] = data.quantity[idx]\n\n            if self.has_stat_profiles:\n                norm_axis = self.stat_scan.geom.axes[\"norm\"]\n                table[\"norm_scan\"] = norm_axis.center.reshape((1, -1))\n                table[\"stat\"] = self.stat.data[idx]\n                table[\"stat_scan\"] = self.stat_scan.data[idx]\n\n            table[\"is_ul\"] = self.is_ul.data[idx]\n            if not self.has_ul:\n                table.remove_columns(\"is_ul\")\n\n        elif format == \"lightcurve\":\n            time_axis = self.geom.axes[\"time\"]\n\n            tables = []\n            for idx, (time_min, time_max) in enumerate(time_axis.iter_by_edges):\n                table_flat = Table()\n                table_flat[\"time_min\"] = [time_min.mjd]\n                table_flat[\"time_max\"] = [time_max.mjd]\n\n                fp = self.slice_by_idx(slices={\"time\": idx})\n                table = fp.to_table(sed_type=sed_type, format=\"gadf-sed\")\n\n                for column in table.columns:\n                    table_flat[column] = table[column][np.newaxis]\n\n                tables.append(table_flat)\n\n            table = vstack(tables)\n\n            # serialize with reference time set to mjd=0.0\n            ref_time = Time(0.0, format=\"mjd\", scale=time_axis.reference_time.scale)\n            table.meta.update(time_ref_to_dict(ref_time, scale=ref_time.scale))\n            table.meta[\"TIMEUNIT\"] = \"d\"\n\n        elif format == \"binned-time-series\":\n            message = (\n                \"Format 'binned-time-series' support a single time axis \"\n                f\"only. Got {self.geom.axes.names}\"\n            )\n\n            if not self.geom.axes.is_unidimensional:\n                raise ValueError(message)\n\n            axis = self.geom.axes.primary_axis\n\n            if not isinstance(axis, TimeMapAxis):\n                raise ValueError(message)\n\n            table = Table()\n            table[\"time_bin_start\"] = axis.time_min\n            table[\"time_bin_size\"] = axis.time_delta\n\n            for quantity in self.all_quantities(sed_type=sed_type):\n                data = getattr(self, quantity, None)\n                if data:\n                    table[quantity] = data.quantity.squeeze()\n        elif format == \"profile\":\n            x_axis = self.geom.axes[\"projected-distance\"]\n\n            tables = []\n            for idx, (x_min, x_max) in enumerate(x_axis.iter_by_edges):\n                table_flat = Table()\n                table_flat[\"x_min\"] = [x_min]\n                table_flat[\"x_max\"] = [x_max]\n                table_flat[\"x_ref\"] = [(x_max + x_min) / 2]\n\n                fp = self.slice_by_idx(slices={\"projected-distance\": idx})\n                table = fp.to_table(sed_type=sed_type, format=\"gadf-sed\")\n\n                for column in table.columns:\n                    table_flat[column] = table[column][np.newaxis]\n\n                tables.append(table_flat)\n\n            table = vstack(tables)\n\n        else:\n            raise ValueError(f\"Not a supported format {format}\")\n\n        if formatted:\n            table = self._format_table(table=table)\n\n        return table\n\n    @staticmethod\n    def _energy_ref_lafferty(model, energy_min, energy_max):\n        \"\"\"Helper for `to_sed_type`.\n\n        Compute energy_ref that the value at energy_ref corresponds\n        to the mean value between energy_min and energy_max.\n        \"\"\"\n        flux = model.integral(energy_min, energy_max)\n        dnde_mean = flux / (energy_max - energy_min)\n        return model.inverse(dnde_mean)\n\n    def _plot_get_flux_err(self, sed_type=None):\n        \"\"\"Compute flux error for given SED type.\"\"\"\n        y_errn, y_errp = None, None\n\n        if \"norm_err\" in self.available_quantities:\n            # symmetric error\n            y_errn = getattr(self, sed_type + \"_err\")",
  "start_line": 450,
  "end_line": 599,
  "total_lines": 872,
  "is_truncated": true
}

{
  "status": "error",
  "error_type": "CommandError",
  "error_message": "",
  "details": {
    "stdout": "",
    "stderr": "",
    "exit_code": 1
  }
}

{
  "status": "ok",
  "stdout": "82:    combined_fp = FluxPoints.from_maps(\n373:        return cls.from_maps(\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "            sed_type = cls._guess_sed_type(table.colnames)\n\n        if sed_type is None:\n            raise ValueError(\"Specifying the SED type is required\")\n\n        if sed_type == \"likelihood\":\n            table = cls._convert_loglike_columns(table)\n            if reference_model is None:\n                reference_model = TemplateSpectralModel(\n                    energy=flat_if_equal(table[\"e_ref\"].quantity),\n                    values=flat_if_equal(table[\"ref_dnde\"].quantity),\n                )\n\n        maps = Maps()\n        table.meta.setdefault(\"SED_TYPE\", sed_type)\n\n        for name in cls.all_quantities(sed_type=sed_type):\n            if name in table.colnames:\n                maps[name] = RegionNDMap.from_table(\n                    table=table, colname=name, format=format\n                )\n\n        meta = cls._get_meta_gadf(table)\n        return cls.from_maps(\n            maps=maps,\n            reference_model=reference_model,\n            meta=meta,\n            sed_type=sed_type,\n            gti=gti,\n        )\n\n    @staticmethod\n    def _get_meta_gadf(table):\n        meta = {}\n        meta.update(table.meta)\n        conf_ul = table.meta.get(\"UL_CONF\")\n        if conf_ul:\n            n_sigma_ul = np.round(stats.norm.isf(0.5 * (1 - conf_ul)), 1)\n            meta[\"n_sigma_ul\"] = n_sigma_ul\n        meta[\"sed_type_init\"] = table.meta.get(\"SED_TYPE\")\n        return meta\n\n    @staticmethod\n    def _format_table(table):\n        \"\"\"Format table.\"\"\"\n        for column in table.colnames:\n            if column.startswith((\"dnde\", \"eflux\", \"flux\", \"e2dnde\", \"ref\")):\n                table[column].format = \".3e\"\n            elif column.startswith(\n                (\"e_min\", \"e_max\", \"e_ref\", \"sqrt_ts\", \"norm\", \"ts\", \"stat\")\n            ):\n                table[column].format = \".3f\"\n\n        return table\n\n    def _guess_format(self):\n        \"\"\"Format of the FluxPoints object.\"\"\"\n        names = self.geom.axes.names\n        if \"time\" in names:\n            return \"lightcurve\"\n        elif \"projected-distance\" in names:\n            return \"profile\"\n        else:\n            return \"gadf-sed\"\n\n    def to_table(self, sed_type=None, format=None, formatted=False):\n        \"\"\"Create table for a given SED type.\n\n        Parameters\n        ----------\n        sed_type : {\"likelihood\", \"dnde\", \"e2dnde\", \"flux\", \"eflux\"}",
  "start_line": 350,
  "end_line": 420,
  "total_lines": 872,
  "is_truncated": false
}

{
  "status": "ok",
  "stdout": "__init__.py\nasmooth.py\ncore.py\nexcess.py\ntests\nts.py\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "stdout": "919:    def from_maps(cls, maps, sed_type=None, reference_model=None, gti=None, meta=None):\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport logging\nfrom copy import deepcopy\nimport numpy as np\nfrom scipy import stats\nfrom scipy.interpolate import interp1d\nfrom scipy.optimize import minimize\nfrom astropy.io import fits\nfrom astropy.io.registry import IORegistryError\nfrom astropy.table import Table, vstack\nfrom astropy.time import Time\nfrom astropy.visualization import quantity_support\nimport matplotlib.pyplot as plt\nfrom gammapy.data import GTI\nfrom gammapy.maps import Map, MapAxis, Maps, RegionNDMap, TimeMapAxis\nfrom gammapy.maps.axes import UNIT_STRING_FORMAT, flat_if_equal\nfrom gammapy.modeling.models import TemplateSpectralModel\nfrom gammapy.modeling.models.spectral import scale_plot_flux\nfrom gammapy.modeling.scipy import stat_profile_ul_scipy\nfrom gammapy.utils.interpolation import interpolate_profile\nfrom gammapy.utils.scripts import make_path\nfrom gammapy.utils.table import table_standardise_units_copy\nfrom gammapy.utils.time import time_ref_to_dict\nfrom ..map.core import DEFAULT_UNIT, FluxMaps\n\n__all__ = [\"FluxPoints\"]\n\nlog = logging.getLogger(__name__)\n\n\ndef squash_fluxpoints(flux_point, axis):\n    \"\"\"Squash a `~FluxPoints` object into one point.\n\n    The log-likelihoods profiles in each bin are summed\n    to compute the resultant quantities. Stat profiles\n    must be present on the `~FluxPoints` object for\n    this method to work.\n    \"\"\"\n    value_scan = flux_point.stat_scan.geom.axes[\"norm\"].center\n    stat_scan = np.sum(flux_point.stat_scan.data, axis=0).ravel()\n    f = interp1d(value_scan, stat_scan, kind=\"quadratic\", bounds_error=False)\n    f = interpolate_profile(value_scan, stat_scan)\n    minimizer = minimize(\n        f,\n        x0=value_scan[int(len(value_scan) / 2)],\n        bounds=[(value_scan[0], value_scan[-1])],\n        method=\"L-BFGS-B\",\n    )\n\n    maps = Maps()",
  "start_line": 1,
  "end_line": 50,
  "total_lines": 872,
  "is_truncated": false
}

{
  "status": "ok",
  "stdout": "40:    stat_scan = np.sum(flux_point.stat_scan.data, axis=0).ravel()\n831:            stat_scan = np.abs(\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/points/core.py",
  "content": "        **kwargs : dict, optional\n            Keyword arguments passed to `~scipy.optimize.brentq`.\n\n        Returns\n        -------\n        flux_points : `~gammapy.estimators.FluxPoints`\n            A new FluxPoints object with modified upper limits.\n\n        Examples\n        --------\n        >>> from gammapy.estimators import FluxPoints\n        >>> filename = '$GAMMAPY_DATA/tests/spectrum/flux_points/binlike.fits'\n        >>> flux_points = FluxPoints.read(filename)\n        >>> flux_points_recomputed = flux_points.recompute_ul(n_sigma_ul=3)\n        >>> print(flux_points.meta[\"n_sigma_ul\"], flux_points.flux_ul.data[0])\n        2.0 [[3.95451985e-09]]\n        >>> print(flux_points_recomputed.meta[\"n_sigma_ul\"], flux_points_recomputed.flux_ul.data[0])\n        3 [[6.22245374e-09]]\n        \"\"\"\n        if not self.has_stat_profiles:\n            raise ValueError(\n                \"Stat profiles not present. Upper limit computation is not possible\"\n            )\n\n        delta_ts = n_sigma_ul**2\n\n        flux_points = deepcopy(self)\n\n        value_scan = self.stat_scan.geom.axes[\"norm\"].center\n        shape_axes = self.stat_scan.geom._shape[slice(3, None)][::-1]\n        for idx in np.ndindex(shape_axes):\n            stat_scan = np.abs(\n                self.stat_scan.data[idx].squeeze() - self.stat.data[idx].squeeze()\n            )\n            flux_points.norm_ul.data[idx] = stat_profile_ul_scipy(\n                value_scan, stat_scan, delta_ts=delta_ts, **kwargs\n            )\n        flux_points.meta[\"n_sigma_ul\"] = n_sigma_ul\n        return flux_points\n\n    def resample_axis(self, axis_new):\n        \"\"\"Rebin the flux point object along the new axis.\n\n        The log-likelihoods profiles in each bin are summed\n        to compute the resultant quantities.\n        Stat profiles must be present on the `~gammapy.estimators.FluxPoints` object for\n        this method to work.\n\n        For now, works only for flat `~gammapy.estimators.FluxPoints`.\n\n        Parameters\n        ----------\n        axis_new : `MapAxis` or `TimeMapAxis`\n            The new axis to resample along\n\n        Returns\n        -------\n        flux_points : `~gammapy.estimators.FluxPoints`\n            A new FluxPoints object with modified axis.\n        \"\"\"\n        if not self.has_stat_profiles:\n            raise ValueError(\"Stat profiles not present, rebinning is not possible\")\n\n        fluxpoints = []\n        for edge_min, edge_max in zip(axis_new.edges_min, axis_new.edges_max):\n            if isinstance(axis_new, TimeMapAxis):\n                edge_min = edge_min + axis_new.reference_time\n                edge_max = edge_max + axis_new.reference_time\n            fp = self.slice_by_coord({axis_new.name: slice(edge_min, edge_max)})\n            fp_new = squash_fluxpoints(fp, axis_new)\n            fluxpoints.append(fp_new)\n\n        return self.__class__.from_stack(fluxpoints, axis=axis_new)",
  "start_line": 800,
  "end_line": 872,
  "total_lines": 872,
  "is_truncated": false
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/map/core.py",
  "content": "            if m.data.dtype is np.dtype(bool):\n                data[name] = m.to_region_nd_map(\n                    region=position, method=\"nearest\", func=np.any\n                )\n            else:\n                data[name] = m.to_region_nd_map(region=position, method=\"nearest\")\n\n        return FluxPoints(\n            data,\n            reference_model=self.reference_model,\n            meta=self.meta.copy(),\n            gti=self.gti,\n        )\n\n    def to_maps(self, sed_type=None):\n        \"\"\"Return maps in a given SED type.\n\n        Parameters\n        ----------\n        sed_type : {\"likelihood\", \"dnde\", \"e2dnde\", \"flux\", \"eflux\"}, optional\n            SED type to convert to. If None, set to `Likelihood`. Default is None.\n\n        Returns\n        -------\n        maps : `Maps`\n            Maps object containing the requested maps.\n        \"\"\"\n        maps = Maps()\n\n        if sed_type is None:\n            sed_type = self.sed_type_init\n\n        for quantity in self.all_quantities(sed_type=sed_type):\n            m = getattr(self, quantity, None)\n            if m is not None:\n                maps[quantity] = m\n\n        return maps\n\n    @classmethod\n    def from_stack(cls, maps, axis, meta=None):\n        \"\"\"Create flux points by stacking list of flux points.\n\n        The first `FluxPoints` object in the list is taken as a reference to infer\n        column names and units for the stacked object.\n\n        Parameters\n        ----------\n        maps : list of `FluxMaps`\n            List of maps to stack.\n        axis : `MapAxis`\n            New axis to create.\n        meta : dict, optional\n            Metadata of the resulting flux points. Default is None.\n\n        Returns\n        -------\n        flux_maps : `FluxMaps`\n            Stacked flux maps along axis.\n        \"\"\"\n        reference = maps[0]\n        data = {}\n\n        for quantity in reference.available_quantities:\n            data[quantity] = Map.from_stack(\n                [_._data[quantity] for _ in maps], axis=axis\n            )\n\n        if meta is None:\n            meta = reference.meta.copy()\n\n        gtis = [_.gti for _ in maps if _.gti]\n\n        if gtis:\n            gti = GTI.from_stack(gtis)\n        else:\n            gti = None\n\n        return cls(\n            data=data, reference_model=reference.reference_model, meta=meta, gti=gti\n        )\n\n    def iter_by_axis(self, axis_name):\n        \"\"\"Create a set of FluxMaps by splitting along an axis.\n\n        Parameters\n        ----------\n        axis_name : str\n             Name of the axis to split on.\n\n        Returns\n        -------\n        flux_maps : `FluxMap`\n            FluxMap iteration.\n\n        \"\"\"\n        split_maps = {}\n        axis = self.geom.axes[axis_name]\n        gti = self.gti\n\n        for amap in self.available_quantities:\n            split_maps[amap] = list(getattr(self, amap).iter_by_axis(axis_name))\n\n        for idx in range(axis.nbin):\n            maps = {}\n            for amap in self.available_quantities:\n                maps[amap] = split_maps[amap][idx]\n                if isinstance(axis, TimeMapAxis):\n                    gti = self.gti.select_time([axis.time_min[idx], axis.time_max[idx]])\n\n            yield self.__class__.from_maps(\n                maps=maps,\n                sed_type=self.sed_type_init,\n                reference_model=self.reference_model,\n                gti=gti,\n                meta=self.meta,\n            )\n\n    @classmethod\n    def from_maps(cls, maps, sed_type=None, reference_model=None, gti=None, meta=None):\n        \"\"\"Create FluxMaps from a dictionary of maps.\n\n        Parameters\n        ----------\n        maps : `Maps`\n            Maps object containing the input maps.\n        sed_type : str, optional\n            SED type of the input maps. If None, set to \"likelihood\". Default is None.\n        reference_model : `~gammapy.modeling.models.SkyModel`, optional\n            Reference model to use for conversions.\n            If None, a model consisting of a point source with a power\n            law spectrum of index 2 is assumed. Default is None.\n        gti : `~gammapy.data.GTI`, optional\n            Maps GTI information. Default is None.\n        meta : `dict`\n            Meta dictionary.\n\n        Returns\n        -------\n        flux_maps : `~gammapy.estimators.FluxMaps`\n            Flux maps object.\n        \"\"\"\n        if sed_type is None:\n            sed_type = cls._guess_sed_type(maps.keys())\n\n        if sed_type is None:\n            raise ValueError(\"Specifying the SED type is required\")\n\n        cls._validate_data(data=maps, sed_type=sed_type)\n",
  "start_line": 800,
  "end_line": 949,
  "total_lines": 1298,
  "is_truncated": true
}

{
  "status": "ok",
  "stdout": "110:class FluxMaps:\n",
  "stderr": "",
  "exit_code": 0
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/map/core.py",
  "content": "class FluxMaps:\n    \"\"\"A flux map / points container.\n\n    It contains a set of `~gammapy.maps.Map` objects that store the estimated\n    flux as a function of energy as well as associated quantities (typically\n    errors, upper limits, delta TS and possibly raw quantities such counts,\n    excesses etc). It also contains a reference model to convert the flux\n    values in different formats. Usually, this should be the model used to\n    produce the flux map.\n\n    The associated map geometry can use a `RegionGeom` to store the equivalent\n    of flux points, or a `WcsGeom`/`HpxGeom` to store an energy dependent flux map.\n\n    The container relies internally on the 'Likelihood' SED type defined in\n    :ref:`gadf:flux-points` and offers convenience properties to convert to\n    other flux formats, namely: ``dnde``, ``flux``, ``eflux`` or ``e2dnde``.\n    The conversion is done according to the reference model spectral shape.\n\n    Parameters\n    ----------\n    data : dict of `~gammapy.maps.Map`\n        The maps dictionary. Expected entries are the following:\n\n        * norm : the norm factor.\n        * norm_err : optional, the error on the norm factor.\n        * norm_errn : optional, the negative error on the norm factor.\n        * norm_errp : optional, the positive error on the norm factor.\n        * norm_ul : optional, the upper limit on the norm factor.\n        * norm_scan : optional, the norm values of the test statistic scan.\n        * stat_scan : optional, the test statistic scan values.\n        * ts : optional, the delta test statistic associated with the flux value.\n        * sqrt_ts : optional, the square root of the test statistic, when relevant.\n        * success : optional, a boolean tagging the validity of the estimation.\n        * n_dof : optional, the number of degrees of freedom used in TS computation\n        * alpha : optional, normalisation factor to accounts for differences between the test region and the background\n        * acceptance_off : optional, acceptance from the off region\n        * acceptance_on : optional, acceptance from the on region\n\n    reference_model : `~gammapy.modeling.models.SkyModel`, optional\n        The reference model to use for conversions. If None, a model consisting\n        of a point source with a power law spectrum of index 2 is assumed.\n    meta : dict, optional\n        Dict of metadata.\n    gti : `~gammapy.data.GTI`, optional\n        Maps GTI information.\n    filter_success_nan : boolean, optional\n        Set fitted norm and error to NaN when the fit has not succeeded.\n    \"\"\"\n\n    _expand_slice = (slice(None), np.newaxis, np.newaxis)\n\n    def __init__(\n        self, data, reference_model, meta=None, gti=None, filter_success_nan=True\n    ):\n        self._data = data\n\n        if isinstance(reference_model, SpectralModel):\n            reference_model = SkyModel(reference_model)\n\n        self._reference_model = reference_model\n\n        if meta is None:\n            meta = {}\n\n        meta.setdefault(\"sed_type_init\", \"likelihood\")\n        self.meta = meta\n        self.gti = gti\n        self._filter_success_nan = filter_success_nan\n\n    @property\n    def filter_success_nan(self):\n        return self._filter_success_nan\n\n    @filter_success_nan.setter\n    def filter_success_nan(self, value):\n        self._filter_success_nan = value\n\n    @property\n    def available_quantities(self):\n        \"\"\"Available quantities.\"\"\"\n        return list(self._data.keys())\n\n    @staticmethod\n    def all_quantities(sed_type):\n        \"\"\"All quantities allowed for a given SED type.\n\n        Parameters\n        ----------\n        sed_type : {\"likelihood\", \"dnde\", \"e2dnde\", \"flux\", \"eflux\"}\n            SED type.\n",
  "start_line": 110,
  "end_line": 200,
  "total_lines": 1298,
  "is_truncated": false
}

{
  "status": "error",
  "error_type": "CommandError",
  "error_message": "",
  "details": {
    "stdout": "",
    "stderr": "",
    "exit_code": 1
  }
}

{
  "status": "error",
  "error_type": "CommandError",
  "error_message": "",
  "details": {
    "stdout": "",
    "stderr": "",
    "exit_code": 1
  },
  "budget_warning": "Only 10 tool call(s) remaining (20/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "stdout": "179:    @property\n187:    @property\n252:    @property\n260:    @property\n265:    @property\n270:    @property\n275:    @property\n280:    @property\n285:    @property\n290:    @property\n311:    @property\n316:    @property\n324:    @property\n329:    @property\n339:    @property\n344:    @property\n349:    @property\n363:    @property\n374:    @property\n386:    @property\n392:    @property\n398:    @property\n429:    @property\n435:    @property\n441:    @property\n448:    @property\n470:    @property\n475:    @property\n480:    @property\n485:    @property\n490:    @property\n495:    @property\n501:    @property\n507:    @property\n513:    @property\n519:    @property\n525:    @property\n531:    @property\n557:    @property\n562:    @property\n568:    @property\n574:    @property\n580:    @property\n586:    @property\n592:    @property\n598:    @property\n604:    @property\n610:    @property\n616:    @property\n622:    @property\n629:    @property\n642:    @property\n655:    @property\n660:    @property\n665:    @property\n670:    @property\n675:    @property\n680:    @property\n685:    @property\n690:    @property\n695:    @property\n700:    @property\n705:    @property\n710:    @property\n715:    @property\n720:    @property\n725:    @property\n730:    @property\n735:    @property\n740:    @property\n745:    @property\n750:    @property\n755:    @property\n",
  "stderr": "",
  "exit_code": 0,
  "budget_warning": "Only 9 tool call(s) remaining (21/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/map/core.py",
  "content": "\n        Returns\n        -------\n        list : list of str\n            All allowed quantities for a given SED type.\n        \"\"\"\n        quantities = []\n        quantities += REQUIRED_MAPS[sed_type]\n        quantities += OPTIONAL_QUANTITIES[sed_type]\n        quantities += OPTIONAL_QUANTITIES_COMMON\n\n        if sed_type == \"likelihood\":\n            quantities += REQUIRED_QUANTITIES_SCAN\n\n        return quantities\n\n    @staticmethod\n    def _validate_data(data, sed_type, check_scan=False):\n        \"\"\"Check that map input is valid and correspond to one of the SED type.\"\"\"\n        try:\n            keys = data.keys()\n            required = set(REQUIRED_MAPS[sed_type])\n        except KeyError:\n            raise ValueError(f\"Unknown SED type: '{sed_type}'\")\n\n        if check_scan:\n            required = required.union(REQUIRED_QUANTITIES_SCAN)\n\n        if not required.issubset(keys):\n            missing = required.difference(keys)\n            raise ValueError(\n                \"Missing data / column for SED type '{}':\" \" {}\".format(\n                    sed_type, missing\n                )\n            )\n\n    # TODO: add support for scan\n    def _check_quantity(self, quantity):\n        if quantity not in self.available_quantities:\n            raise AttributeError(\n                f\"Quantity '{quantity}' is not defined on current flux estimate.\"\n            )\n\n    @staticmethod\n    def _guess_sed_type(quantities):\n        \"\"\"Guess SED type from table content.\"\"\"\n        valid_sed_types = list(REQUIRED_COLUMNS.keys())\n        for sed_type in valid_sed_types:\n            required = set(REQUIRED_COLUMNS[sed_type])\n            if required.issubset(quantities):\n                return sed_type\n\n    @property\n    def is_convertible_to_flux_sed_type(self):\n        \"\"\"Check whether differential SED type is convertible to integral SED type.\"\"\"\n        if self.sed_type_init in [\"dnde\", \"e2dnde\"]:\n            return self.energy_axis.node_type == \"edges\"\n\n        return True\n\n    @property\n    def has_ul(self):\n        \"\"\"Whether the flux estimate has norm_ul defined.\"\"\"\n        return \"norm_ul\" in self._data\n\n    @property\n    def has_any_ts(self):\n        \"\"\"Whether the flux estimate has either sqrt(TS) or test statistic defined.\"\"\"\n        return any(_ in self._data for _ in [\"ts\", \"sqrt_ts\"])\n\n    @property\n    def has_stat_profiles(self):\n        \"\"\"Whether the flux estimate has test statistic profiles.\"\"\"\n        return \"stat_scan\" in self._data\n\n    @property\n    def has_success(self):\n        \"\"\"Whether the flux estimate has the fit status.\"\"\"\n        return \"success\" in self._data\n\n    @property\n    def n_sigma(self):\n        \"\"\"n sigma\"\"\"\n        return self.meta.get(\"n_sigma\", 1)\n\n    @property\n    def n_sigma_ul(self):\n        \"\"\"n sigma UL.\"\"\"\n        return self.meta.get(\"n_sigma_ul\")\n\n    @property\n    def sqrt_ts_threshold_ul(self):\n        \"\"\"sqrt(TS) threshold for upper limits.\"\"\"\n        return self.meta.get(\"sqrt_ts_threshold_ul\", 2)\n\n    @sqrt_ts_threshold_ul.setter\n    def sqrt_ts_threshold_ul(self, value):\n        \"\"\"sqrt(TS) threshold for upper limits.\n\n        Parameters\n        ----------\n        value : int\n            Threshold value in sqrt(TS) for upper limits.\n        \"\"\"\n        self.meta[\"sqrt_ts_threshold_ul\"] = value\n\n        if self.has_any_ts:\n            self.is_ul = self.sqrt_ts < self.sqrt_ts_threshold_ul\n        else:\n            raise ValueError(\"Either TS or sqrt_ts is required to set the threshold\")\n\n    @property\n    def sed_type_init(self):\n        \"\"\"Initial SED type.\"\"\"\n        return self.meta.get(\"sed_type_init\")\n\n    @property\n    def sed_type_plot_default(self):\n        \"\"\"Initial SED type.\"\"\"\n        if self.sed_type_init == \"likelihood\":\n            return \"dnde\"\n\n        return self.sed_type_init\n\n    @property\n    def geom(self):\n        \"\"\"Reference map geometry as a `~gammapy.maps.Geom`.\"\"\"\n        return self.norm.geom\n\n    @property\n    def energy_axis(self):\n        \"\"\"Energy axis as a `~gammapy.maps.MapAxis`.\"\"\"\n        return self.geom.axes[\"energy\"]\n\n    @classproperty\n    def reference_model_default(self):\n        \"\"\"Reference model default as a `~gammapy.modeling.models.SkyModel`.\"\"\"\n        return SkyModel(PowerLawSpectralModel(index=2))\n\n    @property\n    def reference_model(self):\n        \"\"\"Reference model as a `~gammapy.modeling.models.SkyModel`.\"\"\"\n        return self._reference_model\n\n    @property\n    def reference_spectral_model(self):\n        \"\"\"Reference spectral model as a `SpectralModel`.\"\"\"\n        return self.reference_model.spectral_model\n\n    @property",
  "start_line": 200,
  "end_line": 349,
  "total_lines": 1298,
  "is_truncated": true,
  "budget_warning": "Only 8 tool call(s) remaining (22/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "stdout": "558:    def norm(self):\n563:    def norm_err(self):\n569:    def norm_errn(self):\n575:    def norm_errp(self):\n581:    def norm_ul(self):\n587:    def norm_sensitivity(self):\n",
  "stderr": "",
  "exit_code": 0,
  "budget_warning": "Only 7 tool call(s) remaining (23/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/map/core.py",
  "content": "            return self._data[\"sqrt_ts\"]\n        else:\n            with np.errstate(invalid=\"ignore\", divide=\"ignore\"):\n                ts = np.clip(self.ts.data, 0, None)\n                data = np.where(self.norm > 0, np.sqrt(ts), -np.sqrt(ts))\n                return Map.from_geom(geom=self.geom, data=data)\n\n    @property\n    def norm(self):\n        \"\"\"Norm values.\"\"\"\n        return self._filter_convergence_failure(self._data[\"norm\"])\n\n    @property\n    def norm_err(self):\n        \"\"\"Norm error.\"\"\"\n        self._check_quantity(\"norm_err\")\n        return self._filter_convergence_failure(self._data[\"norm_err\"])\n\n    @property\n    def norm_errn(self):\n        \"\"\"Negative norm error.\"\"\"\n        self._check_quantity(\"norm_errn\")\n        return self._data[\"norm_errn\"]\n\n    @property\n    def norm_errp(self):\n        \"\"\"Positive norm error.\"\"\"\n        self._check_quantity(\"norm_errp\")\n        return self._data[\"norm_errp\"]\n\n    @property\n    def norm_ul(self):\n        \"\"\"Norm upper limit.\"\"\"\n        self._check_quantity(\"norm_ul\")\n        return self._data[\"norm_ul\"]\n\n    @property\n    def norm_sensitivity(self):\n        \"\"\"Norm sensitivity.\"\"\"\n        self._check_quantity(\"norm_sensitivity\")\n        return self._data[\"norm_sensitivity\"]\n\n    @property\n    def n_dof(self):\n        \"\"\"Number of degrees of freedom of the fit per energy bin.\"\"\"\n        self._check_quantity(\"n_dof\")\n        return self._data[\"n_dof\"]\n\n    @property\n    def alpha(self):\n        \"\"\"The normalisation, alpha, for differences between the on and off regions.\"\"\"\n        self._check_quantity(\"alpha\")\n        return self._data[\"alpha\"]\n\n    @property\n    def acceptance_on(self):\n        \"\"\"The acceptance in the on region.\"\"\"\n        self._check_quantity(\"acceptance_on\")\n        return self._data[\"acceptance_on\"]\n\n    @property\n    def acceptance_off(self):\n        \"\"\"The acceptance in the off region.\"\"\"\n        self._check_quantity(\"acceptance_off\")\n        return self._data[\"acceptance_off\"]\n\n    @property\n    def dnde_ref(self):\n        \"\"\"Reference differential flux.\"\"\"\n        result = self.reference_spectral_model(self.energy_axis.center)\n        return result[self._expand_slice]\n\n    @property\n    def e2dnde_ref(self):\n        \"\"\"Reference differential flux * energy ** 2.\"\"\"\n        energy = self.energy_axis.center\n        result = self.reference_spectral_model(energy) * energy**2\n        return result[self._expand_slice]\n\n    @property\n    def flux_ref(self):",
  "start_line": 550,
  "end_line": 630,
  "total_lines": 1298,
  "is_truncated": false,
  "budget_warning": "Only 6 tool call(s) remaining (24/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "stdout": "496:    def stat_scan(self):\n",
  "stderr": "",
  "exit_code": 0,
  "budget_warning": "Only 5 tool call(s) remaining (25/30 used). Finalize your edits and call submit_patch soon."
}

{
  "status": "ok",
  "filepath": "gammapy/estimators/map/core.py",
  "content": "    @property\n    def npred_excess_errp(self):\n        \"\"\"Predicted excess counts positive error.\"\"\"\n        return self.npred_excess_ref * self._expand_dims(self.norm_errp.data)\n\n    @property\n    def npred_excess_errn(self):\n        \"\"\"Predicted excess counts negative error.\"\"\"\n        return self.npred_excess_ref * self._expand_dims(self.norm_errn.data)\n\n    @property\n    def npred_excess_ul(self):\n        \"\"\"Predicted excess counts upper limits.\"\"\"\n        return self.npred_excess_ref * self._expand_dims(self.norm_ul.data)\n\n    @property\n    def stat_scan(self):\n        \"\"\"Fit statistic scan value.\"\"\"\n        self._check_quantity(\"stat_scan\")\n        return self._data[\"stat_scan\"]\n\n    @property\n    def dnde_scan_values(self):\n        \"\"\"Fit statistic norm scan values.\"\"\"\n        self._check_quantity(\"dnde_scan_values\")\n        return self._data[\"dnde_scan_values\"]\n\n    @property\n    def stat(self):\n        \"\"\"Fit statistic value.\"\"\"\n        self._check_quantity(\"stat\")",
  "start_line": 480,
  "end_line": 510,
  "total_lines": 1298,
  "is_truncated": false,
  "budget_warning": "Only 4 tool call(s) remaining (26/30 used). Finalize your edits and call submit_patch soon."
}

{
  "error": "Invoking `edit_file()` failed as the following mandatory input parameters are not present:\nold_string\nYou could retry calling this tool, but it is IMPORTANT for you to provide all the mandatory parameters."
}

{
  "error": "Invoking `edit_file()` failed as the following mandatory input parameters are not present:\nold_string\nYou could retry calling this tool, but it is IMPORTANT for you to provide all the mandatory parameters."
}