text:
This is a fixed-text formatted version of a Jupyter notebook

# Maps¶

## Introduction¶

The gammapy.maps submodule contains classes for representing pixilised data on the sky with an arbitrary number of non-spatial dimensions such as energy, time, event class or any possible user-defined dimension (illustrated in the image above). The main Map data structure features a uniform API for WCS as well as HEALPix based images. The API also generalizes simple image based operations such as smoothing, interpolation and reprojection to the arbitrary extra dimensions and makes working with (2 + N)-dimensional hypercubes as easy as working with a simple 2D image. Further information is also provided on the gammapy.maps docs page.

In the following introduction we will learn all the basics of working with WCS based maps. HEALPix based maps will be covered in a future tutorial. We will cover the following topics in order:

Make sure you have worked through the Gammapy overview, because a solid knowledge about working with SkyCoord and Quantity objects as well as Numpy is required for this tutorial.

This notebook is rather lengthy, but getting to know the Map data structure in detail is essential for working with Gammapy and will allow you to fulfill complex analysis tasks with very few and simple code in future!

## Setup¶

[1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import os

[2]:

from astropy import units as u
from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy.convolution import convolve
from gammapy.maps import Map, MapAxis, WcsGeom, WcsNDMap
from gammapy.data import EventList

## Creating WCS Maps¶

### Using Factory Methods¶

Maps are most easily created using the gammapy.maps.Map.create() factory method:

[3]:

m_allsky = Map.create()

Calling Map.create() without any further arguments creates by default an allsky WCS map using a CAR projection, ICRS coordinates and a pixel size of 1 deg. This can be easily checked by printing the .geom attribute of the map:

[4]:

print(m_allsky.geom)

WcsGeom

axes       : ['lon', 'lat']
shape      : (3600, 1800)
ndim       : 2
frame      : icrs
projection : CAR
center     : 0.0 deg, 0.0 deg
width      : 360.0 deg x 180.0 deg
wcs ref    : 0.0 deg, 0.0 deg

The .geom attribute is a Geom object, that defines the basic geometry of the map, such as size of the pixels, width and height of the image, coordinate system etc., but we will learn more about this object later.

Besides the .geom attribute the map has also a .data attribute, which is just a plain numpy.ndarray and stores the data associated with this map:

[5]:

m_allsky.data

[5]:

array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

By default maps are filled with zeros.

The map_type argument can be used to control the pixelization scheme (WCS or HPX).

[6]:

position = SkyCoord(0.0, 5.0, frame="galactic", unit="deg")

# Create a WCS Map
m_wcs = Map.create(binsz=0.1, map_type="wcs", skydir=position, width=10.0)

# Create a HPX Map
m_hpx = Map.create(binsz=0.1, map_type="hpx", skydir=position, width=10.0)

Here is an example that creates a WCS map centered on the Galactic center and now uses Galactic coordinates:

[7]:

skydir = SkyCoord(0, 0, frame="galactic", unit="deg")
m_gc = Map.create(
binsz=0.02, width=(10, 5), skydir=skydir, frame="galactic", proj="TAN"
)
print(m_gc.geom)

WcsGeom

axes       : ['lon', 'lat']
shape      : (500, 250)
ndim       : 2
frame      : galactic
projection : TAN
center     : 0.0 deg, 0.0 deg
width      : 10.0 deg x 5.0 deg
wcs ref    : 0.0 deg, 0.0 deg

In addition we have defined a TAN projection, a pixel size of 0.02 deg and a width of the map of 10 deg x 5 deg. The width argument also takes scalar value instead of a tuple, which is interpreted as both the width and height of the map, so that a quadratic map is created.

### Creating from a Map Geometry¶

As we have seen in the first examples, the Map object couples the data (stored as a numpy.ndarray) with a Geom object. The Geom object can be seen as a generalization of an astropy.wcs.WCS object, providing the information on how the data maps to physical coordinate systems. In some cases e.g. when creating many maps with the same WCS geometry it can be advantegeous to first create the map geometry independent of the map object itsself:

[8]:

wcs_geom = WcsGeom.create(
binsz=0.02, width=(10, 5), skydir=(0, 0), frame="galactic"
)

And then create the map objects from the wcs_geom geometry specification:

[9]:

maps = {}

for name in ["counts", "background"]:
maps[name] = Map.from_geom(wcs_geom)

The Geom object also has a few helpful methods. E.g. we can check whether a given position on the sky is contained in the map geometry:

[10]:

# define the position of the Galactic center and anti-center
positions = SkyCoord([0, 180], [0, 0], frame="galactic", unit="deg")
wcs_geom.contains(positions)

[10]:

array([ True, False])

Or get the image center of the map:

[11]:

wcs_geom.center_skydir

[11]:

<SkyCoord (Galactic): (l, b) in deg
(0., 0.)>

Or we can also retrieve the solid angle per pixel of the map:

[12]:

wcs_geom.solid_angle()

[12]:

$$[[1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~\dots,~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7}],~ [1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~\dots,~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7}],~ [1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~\dots,~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7}],~ \dots,~ [1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~\dots,~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7},~1.2173559 \times 10^{-7}],~ [1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~\dots,~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7},~1.2173376 \times 10^{-7}],~ [1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~\dots,~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7},~1.2173192 \times 10^{-7}]] \; \mathrm{sr}$$

In many analysis scenarios we would like to add extra dimension to the maps to study e.g. energy or time dependency of the data. Those non-spatial dimensions are handled with the MapAxis object. Let us first define an energy axis, with 4 bins:

[13]:

energy_axis = MapAxis.from_bounds(
1, 100, nbin=4, unit="TeV", name="energy", interp="log"
)
print(energy_axis)

MapAxis

name       : energy
unit       : 'TeV'
nbins      : 4
node type  : edges
edges min  : 1.0e+00 TeV
edges max  : 1.0e+02 TeV
interp     : log

Where interp='log' specifies that a logarithmic spacing is used between the bins, equivalent to np.logspace(0, 2, 4). This MapAxis object we can now pass to Map.create() using the axes= argument:

[14]:

m_cube = Map.create(
binsz=0.02, width=(10, 5), frame="galactic", axes=[energy_axis]
)
print(m_cube.geom)

WcsGeom

axes       : ['lon', 'lat', 'energy']
shape      : (500, 250, 4)
ndim       : 3
frame      : galactic
projection : CAR
center     : 0.0 deg, 0.0 deg
width      : 10.0 deg x 5.0 deg
wcs ref    : 0.0 deg, 0.0 deg

Now we see that besides lon and lat the map has an additional axes named energy with 4 bins. The total dimension of the map is now ndim=3.

We can also add further axes by passing a list of MapAxis objects. To demonstrate this we create a time axis with linearly spaced bins and pass both axes to Map.create():

[15]:

time_axis = MapAxis.from_bounds(
0, 24, nbin=24, unit="hour", name="time", interp="lin"
)

m_4d = Map.create(
binsz=0.02, width=(10, 5), frame="galactic", axes=[energy_axis, time_axis]
)
print(m_4d.geom)

WcsGeom

axes       : ['lon', 'lat', 'energy', 'time']
shape      : (500, 250, 4, 24)
ndim       : 4
frame      : galactic
projection : CAR
center     : 0.0 deg, 0.0 deg
width      : 10.0 deg x 5.0 deg
wcs ref    : 0.0 deg, 0.0 deg

The MapAxis object internally stores the coordinates or “position values” associated with every map axis bin or “node”. We distinguish between two node types: edges and center. The node type edges(which is also the default) specifies that the data associated with this axis is integrated between the edges of the bin (e.g. counts data). The node type center specifies that the data is given at the center of the bin (e.g. exposure or differential fluxes).

The edges of the bins can be checked with .edges attribute:

[16]:

energy_axis.edges

[16]:

$$[1,~3.1622777,~10,~31.622777,~100] \; \mathrm{TeV}$$

The numbers are given in the units we specified above, which can be checked again with:

[17]:

energy_axis.unit

[17]:

$$\mathrm{TeV}$$

The centers of the axis bins can be checked with the .center attribute:

[18]:

energy_axis.center

[18]:

$$[1.7782794,~5.6234133,~17.782794,~56.234133] \; \mathrm{TeV}$$

Gammapy Map objects are serialized using the Flexible Image Transport Format (FITS). Depending on the pixelisation scheme (HEALPix or WCS) and presence of non-spatial dimensions the actual convention to write the FITS file is different. By default Gammpy uses a generic convention named gadf, which will support WCS and HEALPix formats as well as an arbitrary number of non-spatial axes. The convention is documented in detail on the Gamma Astro Data Formats page.

Other conventions required by specific software (e.g. the Fermi Science Tools) are supported as well. At the moment those are the following

• fgst-ccube: Fermi counts cube format.

• fgst-ltcube: Fermi livetime cube format.

• fgst-bexpcube: Fermi exposure cube format

• fgst-template: Fermi Galactic diffuse and source template format.

• fgst-srcmap and fgst-srcmap-sparse: Fermi source map and sparse source map format.

The conventions listed above only support an additional energy axis.

Reading FITS files is mainly exposed via the Map.read() method.Let us take a look at a first example:

[19]:

filename = "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz" m_3fhl_gc = Map.read(filename) print(m_3fhl_gc)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : dtype : >i8  If map_type argument is not given when calling read a map object will be instantiated with the pixelization of the input HDU. By default Map.read() will try to find the first valid data hdu in the filename and read the data from there. If mutliple HDUs are present in the FITS file, the desired one can be chosen with the additional hdu= argument: [20]:  m_3fhl_gc = Map.read(filename, hdu="PRIMARY") print(m_3fhl_gc)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : dtype : >i8  In rare cases e.g. when the FITS file is not valid or meta data is missing from the header it can be necessary to modify the header of a certain HDU before creating the Map object. In this case we can use astropy.io.fits directly to read the FITS file: [21]:  filename = ( os.environ["GAMMAPY_DATA"] + "/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz" ) hdulist = fits.open(filename) hdulist.info()  Filename: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 23 (400, 200) float32  And then modify the header keyword and use Map.from_hdulist() to create the Map object after: [22]:  hdulist["PRIMARY"].header["BUNIT"] = "cm2 s" Map.from_hdulist(hdulist=hdulist)  [22]:  WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : cm2 s dtype : >f4  ### Writing Maps¶ Writing FITS files is mainoy exposure via the Map.write() method. Here is a first example: [23]:  m_cube.write("example_cube.fits", overwrite=True)  By default Gammapy does not overwrite files. In this example we set overwrite=True in case the cell gets executed multiple times. Now we can read back the cube from disk using Map.read(): [24]:  m_cube = Map.read("example_cube.fits") print(m_cube)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 250, 4) ndim : 3 unit : dtype : >f4  We can also choose a different FITS convention to write the example cube in a format compatible to the Fermi Galactic diffuse background model: [25]:  m_cube.write("example_cube_fgst.fits", format="fgst-template", overwrite=True)  To understand a little bit better the generic gadf convention we use Map.to_hdulist() to generate a list of FITS HDUs first: [26]:  hdulist = m_4d.to_hdulist(format="gadf") hdulist.info()  Filename: (No file associated with this HDUList) No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 30 (500, 250, 4, 24) float32 1 PRIMARY_BANDS 1 BinTableHDU 33 96R x 7C ['K', 'D', 'D', 'D', 'D', 'D', 'D']  As we can see the HDUList object contains to HDUs. The first one named PRIMARY contains the data array with shape corresponding to our data and the WCS information stored in the header: [27]:  hdulist["PRIMARY"].header  [27]:  SIMPLE = T / conforms to FITS standard BITPIX = -32 / array data type NAXIS = 4 / number of array dimensions NAXIS1 = 500 NAXIS2 = 250 NAXIS3 = 4 NAXIS4 = 24 EXTEND = T WCSAXES = 2 / Number of coordinate axes CRPIX1 = 250.5 / Pixel coordinate of reference point CRPIX2 = 125.5 / Pixel coordinate of reference point CDELT1 = -0.02 / [deg] Coordinate increment at reference point CDELT2 = 0.02 / [deg] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CTYPE1 = 'GLON-CAR' / galactic longitude, plate caree projection CTYPE2 = 'GLAT-CAR' / galactic latitude, plate caree projection CRVAL1 = 0.0 / [deg] Coordinate value at reference point CRVAL2 = 0.0 / [deg] Coordinate value at reference point LONPOLE = 0.0 / [deg] Native longitude of celestial pole LATPOLE = 90.0 / [deg] Native latitude of celestial pole MJDREF = 0.0 / [d] MJD of fiducial time AXCOLS1 = 'E_MIN,E_MAX' INTERP1 = 'log ' AXCOLS2 = 'TIME_MIN,TIME_MAX' INTERP2 = 'lin ' WCSSHAPE= '(500,250,4,24)' BANDSHDU= 'PRIMARY_BANDS' META = '{} ' BUNIT = ''  The second HDU is a BinTableHDU named PRIMARY_BANDS contains the information on the non-spatial axes such as name, order, unit, min, max and center values of the axis bins. We use an astropy.table.Table to show the information: [28]:  Table.read(hdulist["PRIMARY_BANDS"])  [28]:  Table length=96 CHANNELENERGYE_MINE_MAXTIMETIME_MINTIME_MAX TeVTeVTeVhhh int64float64float64float64float64float64float64 01.7782794100389231.03.16227766016837950.50.01.0 15.6234132519034923.162277660168379510.0000000000000020.50.01.0 217.78279410038923610.00000000000000231.6227766016838030.50.01.0 356.23413251903491531.622776601683803100.000000000000040.50.01.0 41.7782794100389231.03.16227766016837951.51.02.0 55.6234132519034923.162277660168379510.0000000000000021.51.02.0 617.78279410038923610.00000000000000231.6227766016838031.51.02.0 756.23413251903491531.622776601683803100.000000000000041.51.02.0 81.7782794100389231.03.16227766016837952.52.03.0 ..................... 8617.78279410038923610.00000000000000231.62277660168380321.521.022.0 8756.23413251903491531.622776601683803100.0000000000000421.521.022.0 881.7782794100389231.03.162277660168379522.522.023.0 895.6234132519034923.162277660168379510.00000000000000222.522.023.0 9017.78279410038923610.00000000000000231.62277660168380322.522.023.0 9156.23413251903491531.622776601683803100.0000000000000422.522.023.0 921.7782794100389231.03.162277660168379523.523.024.0 935.6234132519034923.162277660168379510.00000000000000223.523.024.0 9417.78279410038923610.00000000000000231.62277660168380323.523.024.0 9556.23413251903491531.622776601683803100.0000000000000423.523.024.0 Maps can be serialized to a sparse data format by calling write with sparse=True. This will write all non-zero pixels in the map to a data table appropriate to the pixelization scheme. [29]:  m = Map.create(binsz=0.1, map_type="wcs", width=10.0) m.write("file.fits", hdu="IMAGE", sparse=True, overwrite=True) m = Map.read("file.fits", hdu="IMAGE", map_type="wcs")  ## Accessing Data¶ ### How to get data values¶ All map objects have a set of accessor methods, which can be used to access or update the contents of the map irrespective of its underlying representation. Those accessor methods accept as their first argument a coordinate tuple containing scalars, list, or numpy.ndarray with one tuple element for each dimension. Some methods additionally accept a dict or MapCoord argument, of which both allow to assign coordinates by axis name. Let us first begin with the .get_by_idx() method, that accepts a tuple of indices. The order of the indices corresponds to the axis order of the map: [30]:  m_gc.get_by_idx((50, 30))  [30]:  array([0.], dtype=float32)  Important: Gammapy uses a reversed index order in the map API with the longitude axes first. To achieve the same by directly indexing into the numpy array we have to call: [31]:  m_gc.data[([30], [50])]  [31]:  array([0.], dtype=float32)  To check the order of the axes you can always print the .geom attribute: [32]:  print(m_gc.geom)  WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 frame : galactic projection : TAN center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg  To access values directly by sky coordinates we can use the .get_by_coord() method. This time we pass in a dict, specifying the axes names corresponding to the given coordinates: [33]:  m_gc.get_by_coord({"lon": [0, 180], "lat": [0, 0]})  [33]:  array([ 0., nan])  The units of the coordinates are assumed to be in degrees in the coordinate system used by the map. If the coordinates do not correspond to the exact pixel center, the value of the nearest pixel center will be returned. For positions outside the map geometry np.nan is returned. The coordinate or idx arrays follow normal Numpy broadcasting rules. So the following works as expected: [34]:  lons = np.linspace(-4, 4, 10) m_gc.get_by_coord({"lon": lons, "lat": 0})  [34]:  array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)  Or as an even more advanced example, we can provide lats as column vector and broadcasting to a 2D result array will be applied: [35]:  lons = np.linspace(-4, 4, 8) lats = np.linspace(-4, 4, 8).reshape(-1, 1) m_gc.get_by_coord({"lon": lons, "lat": lats})  [35]:  array([[nan, nan, nan, nan, nan, nan, nan, nan], [nan, nan, nan, nan, nan, nan, nan, nan], [ 0., 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0., 0.], [nan, nan, nan, nan, nan, nan, nan, nan], [nan, nan, nan, nan, nan, nan, nan, nan]])  ### Indexing and Slicing Sub-Maps¶ When you have worked with Numpy arrays in the past you are probably familiar with the concept of indexing and slicing into data arrays. To support slicing of non-spatial axes of Map objects, the Map object has a .slice_by_idx() method, which allows to extract sub-maps from a larger map. The following example demonstrates how to get the map at the energy bin number 3: [36]:  m_sub = m_cube.slice_by_idx({"energy": 3}) print(m_sub)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 unit : dtype : >f4  Note that the returned object is again a Map with updated axes information. In this case, because we extracted only a single image, the energy axes is dropped from the map. To extract a sub-cube with a sliced energy axes we can use a normal slice() object: [37]:  m_sub = m_cube.slice_by_idx({"energy": slice(1, 3)}) print(m_sub)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 250, 2) ndim : 3 unit : dtype : >f4  Note that the returned object is also a Map object, but this time with updated energy axis specification. Slicing of multiple dimensions is supported by adding further entries to the dict passed to .slice_by_idx() [38]:  m_sub = m_4d.slice_by_idx({"energy": slice(1, 3), "time": slice(4, 10)}) print(m_sub)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy', 'time'] shape : (500, 250, 2, 6) ndim : 4 unit : dtype : float32  For convenience there is also a .get_image_by_coord() method which allows to access image planes at given non-spatial physical coordinates. This method also supports Quantity objects: [39]:  image = m_4d.get_image_by_coord({"energy": 4 * u.TeV, "time": 5 * u.h}) print(image.geom)  WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg  ### Iterating by image¶ For maps with non-spatial dimensions the ~Map.iter_by_image method can be used to loop over image slices. The image plane index idx is returned in data order, so that the data array can be indexed directly. Here is an example for an in-place convolution of an image using astropy.convolution.convolve to interpolate NaN values: [40]:  axis1 = MapAxis([1, 10, 100], interp="log", name="energy") axis2 = MapAxis([1, 2, 3], interp="lin", name="time") m = Map.create(width=(5, 3), axes=[axis1, axis2], binsz=0.1) m.data[:, :, 15:18, 20:25] = np.nan for img, idx in m.iter_by_image(): kernel = np.ones((5, 5)) m.data[idx] = convolve(img, kernel) assert not np.isnan(m.data).any()  ## Modifying Data¶ ### How to set data values¶ To modify and set map data values the Map object features as well a .set_by_idx() method: [41]:  m_cube.set_by_idx(idx=(10, 20, 3), vals=42)  here we check that data have been updated: [42]:  m_cube.get_by_idx((10, 20, 3))  [42]:  array([42.], dtype=float32)  Of course there is also a .set_by_coord() method, which allows to set map data values in physical coordinates. [43]:  m_cube.set_by_coord({"lon": 0, "lat": 0, "energy": 2 * u.TeV}, vals=42)  Again the lon and lat values are assumed to be given in degrees in the coordinate system used by the map. For the energy axis, the unit is the one specified on the axis (use m_cube.geom.axes[0].unit to check if needed…). All .xxx_by_coord() methods accept SkyCoord objects as well. In this case we have to use the skycoord keyword instead of lon and lat: [44]:  skycoords = SkyCoord([1.2, 3.4], [-0.5, 1.1], frame="galactic", unit="deg") m_cube.set_by_coord({"skycoord": skycoords, "energy": 2 * u.TeV}, vals=42)  ### Filling maps from event lists¶ This example shows how to fill a counts cube from an event list: [45]:  energy_axis = MapAxis.from_bounds( 10.0, 2e3, 12, interp="log", name="energy", unit="GeV" ) counts_3d = WcsNDMap.create( binsz=0.1, width=10.0, skydir=(0, 0), frame="galactic", axes=[energy_axis] ) events = EventList.read( "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.fits.gz"
)

counts_3d.write("ccube.fits", format="fgst-ccube", overwrite=True)

Alternatively you can use the fill_events method:

[46]:

counts_3d = WcsNDMap.create(
binsz=0.1, width=10.0, skydir=(0, 0), frame="galactic", axes=[energy_axis]
)

counts_3d.fill_events(events)

If you have a given map already, and want to make a counts image with the same geometry (not using the pixel data from the original map), you can also use the fill_events method.

[47]:

events = EventList.read(
"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.fits.gz" ) reference_map = Map.read( "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz"
)
counts = Map.from_geom(reference_map.geom)
counts.fill_events(events)

It works for IACT and Fermi-LAT events, for WCS or HEALPix map geometries, and also for extra axes. Especially energy axes are automatically handled correctly.

### Filling maps from interpolation¶

Maps support interpolation via the ~Map.interp_by_coord and ~Map.interp_by_pix methods. Currently the following interpolation methods are supported:

• nearest : Return value of nearest pixel (no interpolation).

• linear : Interpolation with first order polynomial. This is the only interpolation method that is supported for all map types.

• quadratic : Interpolation with second order polynomial.

• cubic : Interpolation with third order polynomial.

Note that quadratic and cubic interpolation are currently only supported for WCS-based maps with regular geometry (e.g. 2D or ND with the same geometry in every image plane). linear and higher order interpolation by pixel coordinates is only supported for WCS-based maps.

In the following example we create a new map and fill it by interpolating another map:

[48]:

# read map
filename = "$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz" m_iem_gc = Map.read(filename) # create new geometry skydir = SkyCoord(266.4, -28.9, frame="icrs", unit="deg") wcs_geom_cel = WcsGeom.create( skydir=skydir, binsz=0.1, frame="icrs", width=(8, 4) ) # create new empty map from geometry m_iem_10GeV = Map.from_geom(wcs_geom_cel) coords = m_iem_10GeV.geom.get_coord() # fill new map using interpolation m_iem_10GeV.data = m_iem_gc.interp_by_coord( {"skycoord": coords.skycoord, "energy_true": 10 * u.GeV}, method="linear", fill_value=np.nan, )  ### Interpolating onto a different geometry¶ For 3d geometries this operation can be performed directly using the .interp_to_geom() method. This is very useful, ex: while using map arithmetics. [49]:  # create new geometry energy_axis = MapAxis.from_bounds( 10.0, 2e3, 6, interp="log", name="energy_true", unit="GeV" ) skydir = SkyCoord(266.4, -28.9, frame="icrs", unit="deg") wcs_geom_3d = WcsGeom.create( skydir=skydir, binsz=0.1, frame="icrs", width=(8, 4), axes=[energy_axis] ) # create the interpolated map m_iem_interp = m_iem_gc.interp_to_geom( wcs_geom_3d, preserve_counts=False, method="linear", fill_value=np.nan ) print(m_iem_interp)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy_true'] shape : (80, 40, 6) ndim : 3 unit : 1 / (cm2 MeV s sr) dtype : float64  Note that preserve_counts option should be true if the map is an integral quantity (e.g. counts) and false if the map is a differential quantity (e.g. intensity). ## Maps operations¶ ### Basic operators¶ One can perform simple arithmetic on maps using the +, -, *, / operators, this works only for maps with the same geometry: [50]:  iem_plus_iem = m_iem_10GeV + m_iem_10GeV iem_minus_iem = m_iem_10GeV - m_iem_10GeV  These operations can be applied between a Map and a scalar in that specific order: [51]:  iem_times_two = m_iem_10GeV * 2 # iem_times_two = 2 * m_iem_10GeV # this won't work  The logic operators can also by applied on maps (the result is a map of boolean type): [52]:  is_null = iem_minus_iem == 0 print(is_null)  WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (80, 40) ndim : 2 unit : dtype : bool  Here we check that the result is True for all the well-defiend pixels (not NaN): [53]:  np.all(is_null.data[~np.isnan(iem_minus_iem)])  [53]:  True  ### Cutouts¶ The WCSNDMap objects features a .cutout() method, which allows you to cut out a smaller part of a larger map. This can be useful, e.g. when working with allsky diffuse maps. Here is an example: [54]:  position = SkyCoord(0, 0, frame="galactic", unit="deg") m_iem_cutout = m_iem_gc.cutout(position=position, width=(4 * u.deg, 2 * u.deg))  The returned object is again a Map object with udpated WCS information and data size. As one can see the cutout is automatically applied to all the non-spatial axes as well. The cutout width is given in the order of (lon, lat) and can be specified with units that will be handled correctly. ## Visualizing and Plotting¶ All map objects provide a plot method for generating a visualization of a map. This method returns figure, axes, and image objects that can be used to further tweak/customize the image. The plot method should be used with 2D maps, while 3D maps can be displyed with the .plot_interative() or .plot_grid() methods. ### Image Plotting¶ For debugging and inspecting the map data it is useful to plot ot visualize the images planes contained in the map. [55]:  filename = "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz"

After reading the map we can now plot it on the screen by calling the .plot() method:

[56]:

m_3fhl_gc.plot();

We can easily improve the plot by calling Map.smooth() first and providing additional arguments to .plot(). Most of them are passed further to plt.imshow():

[57]:

smoothed = m_3fhl_gc.smooth(width=0.2 * u.deg, kernel="gauss")

/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
return super().imshow(X, *args, origin=origin, **kwargs)

We can use the plt.rc_context() context manager to further tweak the plot by adapting the figure and font size:

[58]:

rc_params = {"figure.figsize": (12, 5.4), "font.size": 12}
with plt.rc_context(rc=rc_params):
smoothed = m_3fhl_gc.smooth(width=0.2 * u.deg, kernel="gauss")

/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
return super().imshow(X, *args, origin=origin, **kwargs)

### Cube plotting¶

For maps with non-spatial dimensions the Map object features an interactive plotting method, that works in jupyter notebooks only (Note: it requires the package ipywidgets to be installed). We first read a small example cutout from the Fermi Galactic diffuse model and display the data cube by calling .plot_interactive():

[59]:

rc_params = {
"figure.figsize": (12, 5.4),
"font.size": 12,
"axes.formatter.limits": (2, -2),
}

Now you can use the interactive slider to select an energy range and the corresponding image is diplayed on the screen. You can also use the radio buttons to select your preferred image stretching. We have passed additional keywords using the rc_params argument to improve the figure and font size. Those keywords are directly passed to the plt.rc_context() context manager.

Additionnaly all the slices of a 3D Map can be displayed using the .plot_grid() method. By default the colorbars bounds of the subplots are not the same, we can make them consistant using the vmin and vmax options:

[60]:

counts_3d.plot_grid(ncols=4, figsize=(16, 12), vmin=0, vmax=100, stretch="log");

/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
return super().imshow(X, *args, origin=origin, **kwargs)

[ ]:


text:
# Twelve Years a Slave (2013 film) Summary and Analysis of Part 2

Summary

When they arrive at their destination, Clemens' master is waiting for him and demands that he return to the plantation. Clemens runs to his master, pretending to be much more grateful and less intelligent than he has shown himself to be to Solomon, leaving Solomon behind. We see a flashback of Solomon going into a store with his family in Saratoga, where Solomon is buying a carryall for his wife's travels. While they are shopping, a black man comes into the store, followed by his master, a white man, who scolds him for wandering off.

Back in the present, a slaver named Freeman calls the remaining slaves' names and slaps Solomon when he does not immediately answer to his new slave name, "Platt." At Freeman's offices, Freeman tries to sell his slaves to some white buyers. When a man named Ford is interested in Solomon (for $1000) and Eliza (for$700), Eliza begs Freeman to sell her children with her, but he swiftly sells her son (for \$600) to a different buyer. Ford tries to buy Eliza's daughter, but Freeman insists that she is more valuable, because she is half-white, and that he cannot sell her with her mother. Eliza is separated from both her children. When she leaves the room, screaming, Solomon begins to play the violin.

Ford brings Eliza and Solomon to his home, where his wife looks at the grieving Eliza. "Something to eat and some rest. Your children will soon be forgotten," she says, going into the house.

We see slave men preparing to work, and a white slave-handler and carpenter named John Tibeats introduces himself and Mr. Chapin, an overseer on the plantation. He tells the slave that they must refer to both of them as "master" and invites them all to clap their hands while he sings a racist song about killing black men. The slaves cut down trees for their work that day. We see Ford delivering a Christian sermon to his slaves.

One day, Solomon and the other slave men run into a group of Native Americans, with whom they spend a brief moment of peace. Solomon watches one of the men play a stringed instrument and seems to think of his violin.

Solomon tells Tibeats and Ford about an idea of how to transport lumber more quickly on the river. "Are you an engineer or a nigger?" Tibeats asks, patronizingly. Ford, on the other hand, wants to hear him out, and is impressed as Solomon tells them he once hired his own group to transport items from Lake Champlain to Troy, New York. He tells Solomon to collect a team and put his plan into action.

Solomon puts his plan into motion and succeeds, much to Tibeats' chagrin. Ford is impressed and rewards him with a violin. "I hope it brings us both much joy over the years," Ford says, as Solomon goes back to the slave quarters. There, Eliza is sobbing, and Solomon scolds her, saying, "You let yourself be overcome by sorrow, you'll drown in it." She asks him whether he mourns his children, and he says that he misses them a great deal.

"Master Ford is a decent man," Solomon protests, but Eliza insists, "He is a slaver." Solomon insists that he does what he needs to do to survive, but Eliza tells him, "You are no better than prized livestock." The scene shifts to Ford delivering a sermon. While Eliza sobs, Mrs. Ford whispers, "I cannot have that kind of depression about."

Later, Tibeats approaches Solomon and insults him. "You are a dog," Tibeats says, barking at him and telling him to work harder on the house he is erecting. Suddenly, Solomon sees Eliza being dragged off. We see a flashback of Eliza telling a story of her old life on a plantation where she was treated well and even bore a daughter with her master. When that master got older and lost power in the household, she says, the master's resentful daughter brought her and her children to the city under the pretense that they were going to be freed.

Back in the present, Solomon works on the clapboards on the house, when he is once again visited by Tibeats. When Solomon speaks back to Tibeats, Tibeats tells him to strip, but Solomon refuses. The two men fight, with Solomon eventually beating Tibeats. Chapin approaches and tells Solomon that if he runs, he cannot protect him, implying that he will go to Ford to sort out the conflict.

Tibeats hangs Solomon from a tree, but Chapin approaches and warns him that he has no claim to Solomon's life, as Ford is his owner.

Analysis

Solomon's fate seems only to get bleaker and bleaker as he descends deeper and deeper into the world of slavery. After getting cut off from his only friend, Clemens, when the latter returns to his master, Solomon is left alone in the care of a particularly ruthless slaver named Freeman, who forces Solomon to take a new name. Solomon is brought back to Freeman's office where he is sold and witnesses the tragic separation of Eliza from her children. The film plunges the viewer into the inhumanity and horrors of slavery.

One of the more horrific elements of the narrative is the fact that while some white characters are more sympathetic towards the plight of the slaves, they do nothing to stop the system itself. Presented alongside a character like the bigoted and evil Freeman, Mr. Ford and his wife seem to have some goodness in them, in that they try to pity the people they have enslaved. When Mrs. Ford sees Eliza crying, she appears to feel sorry for the woman's loss of her children, but soon follows up her sentiment of sympathy with the statement, "Your children will soon be forgotten." In this, we see that Mrs. Ford's cruelty is as damaging and unconscious as the ruthless Freeman, if more sympathetically packaged. While there seems to be a spectrum of cruelty, all of the white characters in the film are presented as complicit in upholding the institution of slavery.

McQueen juxtaposes the horrors of the human slave trade with the beauty of the natural world throughout the film. On the boat, he films the water of the river rushing through the machinery propelling the boat forward. On the plantation, he shoots beautifully lit shots of the majestic and mysterious trees on the property. In these images, McQueen juxtaposes images of the natural world with images of human labor. We see the mechanics of the boat moving towards its destination, all part of an industry dedicated to the subjugation of slaves. Immediately after showing the trees on the Ford property, McQueen shows a line of slave men preparing to work. While the message of these images is implicit, they seem to suggest that slavery and colonialism both concern the subjugation of nature and human life; McQueen reminds us that the enjoyment of the natural world within American colonial structures comes at a high human price.

Solomon finds himself at odds with some of the other slaves, when it comes to his desire to win favor with his master. After convincing Ford to let him pursue an engineering project that impresses his master, Solomon scolds Eliza, his fellow slave, for mourning the loss of her children. In this moment, she reminds him that, while Ford may be showing kindness to him, he is a slaver, and he thinks of Solomon as little more than "prized livestock." Solomon, a Northerner who is used to doing what he can to climb the ladder of white society, finds himself brought back down to earth in this moment, reminded that even if the white people on the plantation seem to believe in him, they view him with a racist condescension.

The tragedy is that neither Solomon nor Eliza's tactics are effective in resisting the indignity and tragedy of slavery. By openly expressing her grief at having lost her children, Eliza makes herself vulnerable to the whims of her master. Similarly, by keeping his mouth shut and simply trying to excel, Solomon becomes a target for Tibeats, who views Solomon as a threat to his own dignity and intellect. Solomon and Eliza's clash is all the more tragic because neither of them is able to maintain any sense of dignity under slavery, because slavery does not view them as being dignified human beings.

format: plain


text:
# Real-time Likelihood-free Inference of Roman Binary Microlensing Events with Amortized Neural Posterior Estimation

### Abstract

Fast and automated inference of binary-lens, single-source (2L1S) microlensing events with sampling-based Bayesian algorithms (e.g., Markov Chain Monte Carlo, MCMC) is challenged on two fronts: the high computational cost of likelihood evaluations with microlensing simulation codes, and a pathological parameter space where the negative-log-likelihood surface can contain a multitude of local minima that are narrow and deep. Analysis of 2L1S events usually involves grid searches over some parameters to locate approximate solutions as a prerequisite to posterior sampling, an expensive process that often requires human-in-the- loop domain expertise. As the next-generation, space-based microlensing survey with the Roman Space Telescope is expected to yield thousands of binary microlensing events, a new fast and automated method is desirable. Here, we present a likelihood- free inference approach named amortized neural posterior estimation, where a neural density estimator (NDE) learns a surrogate posterior $hatp(boldsymbolþeta | boldsymbolx)$ as an observation-parameterized conditional probability distribution, from pre-computed simulations over the full prior space. Trained on 291,012 simulated Roman-like 2L1S simulations, the NDE produces accurate and precise posteriors within seconds for any observation within the prior support without requiring a domain expert in the loop, thus allowing for real-time and automated inference. We show that the NDE also captures expected posterior degeneracies. The NDE posterior could then be refined into the exact posterior with a downstream MCMC sampler with minimal burn-in steps.

Publication
aj


text:
Publications: Abstract

# General Adversaries in Unconditional Multi-Party Computation

## Matthias Fitzi and Martin Hirt and Ueli Maurer

We consider a generalized adversary model for unconditionally secure multi-party computation. The adversary can actively corrupt (i.e. take full control over) a subset \$\ClassA\subseteq\PS\$ of the players, and, additionally}, can passively corrupt (i.e. read the entire information of) another subset \$\ClassP\subseteq\PS\$ of the players. The adversary is characterized by a generalized adversary structure, i.e. a set of pairs \$(\ClassA,\ClassP)\$, where he may select one arbitrary pair from the structure and corrupt the players accordingly. This generalizes the classical threshold results of Ben-Or, Goldwasser and Wigderson, Chaum, Crépeau, and Damg{å}rd, and Rabin and Ben-Or, and the non-threshold results of Hirt and Maurer.

The generalizations and improvements on the results of Hirt and Maurer are three-fold: First, we generalize their model by considering mixed (active and passive) non-threshold adversaries and characterize completely the adversary structures for which unconditionally secure multi-party computation is possible, for four different models: Perfect security with and without broadcast, and unconditional security (with negligible error probability) with and without broadcast. All bounds are tight. Second, some of their protocols have complexity super-polynomial in the size of the adversary structure; we reduce the complexity to polynomial. Third, we prove the existence of adversary structures for which no polynomial (in the number of players) protocols exist.

The following two implications illustrate the usefulness of these results: The most powerful adversary that is unconditionally tolerated by previous protocols among three players is the one that passively corrupts one arbitrary player; using our protocols one can unconditionally tolerate an adversary that either passively corrupts the first player, or actively\/} corrupts the second or the third player.

Moreover, in a setting with arbitrarily many cheating players who want to compute an agreed function with the help of a trusted party, we can relax the trust requirement into this helping party: Without support from the cheating players the helping party obtains no information about the honest players' inputs and outputs.

Keywords: General adversaries, mixed model, multi-party computation, unconditional security.


text:
# Useful identities from Euler’s formula

Euler’s formula:

$e^{i\theta} =\cos\theta + i\sin\theta$

From the Euler’s formula, we can derive that

$e^{ik\pi} = \begin{cases} 1 & \text{if } k \text{ is even}\\ -1 & \text{if } k \text{ is odd} \end{cases}$

If we write the above equation in reverse, we get

$-1=e^{i(2k+1)\pi}$

$1=e^{i2k\pi}$

These identities are useful for quick manipulations of complex numbers. Below are two examples:

Example 1: Solve $z^n = 1$.

Solution

$z = 1^{1/n} = e^{i\frac{2k\pi}{n}} \qquad k = 0,1,...,n-1$

Example 2: Solve $z^3=-8$.

Solution

$z = (-8)^{1/3} = 2(-1)^{1/3} = 2e^{i\frac{(2k+1)\pi}{3}} \qquad k =0, 1, 2$

• $k=0 \Rightarrow z=1+\sqrt{3}i$
• $k=1 \Rightarrow z=-2$
• $k=2 \Rightarrow z=1-\sqrt{3}i$

Note that we can also write $k=-1, 0, 1$. On the trigonometric circle, $k = 2$ and $k = -1$ yield the same angle.


text:
# Ways of merging two incomparable sorted lists of elements keeping their relative ordering

Suppose that, for a real application, I have ended up with a sorted list A = {$a_1, a_2, ..., a_{|A|}$} of elements of a certain kind (say, Type-A), and another sorted list B = {$b_1, b_2, ..., b_{|B|}$} of elements of a different kind (Type-B), such that Type-A elements are only comparable with Type-A elements, and likewise for Type-B.

At this point I seek to count the following: in how many ways can I merge both lists together, in such a way that the relative ordering of Type-A and Type-B elements, respectively, is preserved? (i.e. that if $P_M(x)$ represents the position of an element of A or B in the merged list, then $P_M(a_i)<P_M(a_j)$ and $P_M(b_i)<P_M(b_j)$ for all $i<j$)

I've tried to figure this out constructively by starting with an empty merged list and inserting elements of A or B one at a time, counting in how many ways each insertion can be done, but since this depends on the placement of previous elements of the same type, I've had little luck so far. I also tried explicitly counting all possibilities for different (small) lengths of A and B, but I've been unable to extract any potential general principle in this way.

The merged list has length $|A|+|B|$. Once you know which $|A|$ positions in it are occupied by the elements of list A, you also know exactly how the whole thing has to be ordered, since the internal orders of the elements of A and the elements of B are already known. Thus, there are $$\binom{|A|+|B|}{|A|}=\binom{|A|+|B|}{|B|}$$ possible merged lists, one for each choice of $|A|$ positions for the elements of list A.

• Great answer! I've made an impl in C++ that generates all possible interleaved combinations of two lists, while maintaining relative order for anyone who is interested: github.com/domfarolino/algorithms/tree/master/src/… Oct 22, 2018 at 20:37
• Can this further generalize to 3 lists? To 4? To N? May 27, 2021 at 15:54
• Indeed, I believe this can further generalize to N lists, by viewing the overall operation as a series of (N-1) merges, so the total number of merged lists is computed as the product of (N-1) Binomial Coeffecients, where the i-th expression computes the selection of the positions of the i-th lists's elements from the available positions left after the first (i-1) lists have already been merged. May 27, 2021 at 21:38
• @BryanPendleton: Yes: it’s a multinomial coefficient, $$\binom{|A_1|+|A_2|+\ldots+|A_N|}{A_1,A_2,\ldots,A_N}\,,$$ which can be expanded as a product of binomial coefficients as shown at the link. May 28, 2021 at 3:57
• Great answer. Since we have to keep the relative ordering. This question is equal to "How many different ways are there to select n positions from m + n positions?“ Sep 17, 2021 at 5:47

$\frac{(m+n)!}{(m)!(n)!}$; where m,n is the length of the arrays. $(m+n)!$ is the total number of ways of arranging without any restriction. The condition given that you can not change the order of each array so you should not change the order of each array , so we need to divide all the permutations of each array!

For merging N sorted lists, here is a good way to see that the solution is $$\frac{(|A_1|+\dots|+|A_N|)!}{|A_1|!\dots |A_N|!}$$

All the $$(|A_1|+\dots|+|A_N|)$$ elements can be permuted in $$(|A_1|+\dots|+|A_N|)!$$ ways. Among these, any solution which has the ordering of the $$A_1$$ elements different from the given order have to be thrown out. If you keep the positions of everyone except $$A_1$$ fixed, you can generate $$|A_1|!$$ permutations by shuffling $$A_1$$. Only 1 among these $$|A_1|!$$ solutions is valid. Hence the number solutions containing the right order of $$A_1$$ is $$\frac{(|A_1|+\dots|+|A_N|)!}{|A_1|!}$$

Now repeat the argument for $$A_2$$, $$A_2, \dots A_N$$.


text:
# ADM decomposition of the general scalar tensor theory Lagrangian

I have question about ADM decomposition of some general scalar-tensor theory of gravity.

Starting with ADM form of the metric:

$$ds^2=-N^2dt^2+h_{ij}(dx^i+N^idt)(dx^j+N^jdt)$$

provided with extrinsic curvature:

$$K_{ij}=\frac{1}{2N}(\dot{h_{ij}}-D_iN_j-D_jN_i)$$, where $$N$$ is lapse function and $$N^i$$'s are shift vector components. Timelike normal vector to the hypersurface is denoted as $$n_a$$.

I'm considering given action:

$$\int d^4x\sqrt{-g}[f(\phi) R-\nabla_\mu\phi\nabla^\mu\phi + U(\phi)]$$

and I want to recast this Lagrangian in the 3+1 form which is suitable for discussing Hamiltonian formulation of this theory.

Kinetic term $$(X=\nabla_\mu\phi\nabla^\mu\phi)$$ decomposes as: $$X=-A_*^2+D^i\phi D_i\phi$$

where: $$A_*=n^\mu\nabla_\mu=\frac {1}{N}(\dot{\phi}-N^iD_i\phi)$$ and $$D^i$$ is 3d covariant derivative associated with metric $$h_{ab}$$ ($$h^a_b\nabla_a=D_b$$).

However, i have trouble with the first part of the action:

$$\int d^4x\sqrt{-g}f(\phi) R=\int dt\int d^3x \sqrt{-h}N\Big[f(\phi)\big(R^{(3)}+K_{ij}K^{ij}-K^2+2\nabla_{\mu}(n^\mu\nabla_\nu n^\nu-n^\nu\nabla_\nu n^\mu)\big)\Big]$$.

First term is left as it is, while term involving derivatives of normal vector $$n^\mu$$ needs to be integrated by parts (dropping boundary term) - normally in GR this is total divergence and is discarded (this is not the case here):

$$\int dt \int d^3x\sqrt{-h}N f(\phi) 2\nabla_{\mu}(n^\mu\nabla_\nu n^\nu-n^\nu\nabla_\nu n^\mu) \\=-2\int dt \int d^3x\sqrt{-h}N(n^\mu\nabla_\nu n^\nu-n^\nu\nabla_\nu n^\mu)\nabla_\mu f(\phi) \\= -2\int dt \int d^3x\sqrt{-h}N(n^\mu\nabla_\nu n^\nu)f_\phi\nabla_\mu\phi+2\int dt \int d^3x\sqrt{-h}N(n^\nu\nabla_\nu n^\mu)f_\phi \nabla_\mu \phi\\=-2\int dt \int d^3x\sqrt{-h}N(n^\mu Kf_\phi\nabla_\mu\phi-n^\nu(\nabla_\nu n^\mu) f_\phi \nabla_\mu \phi)=-2\int dt \int d^3x\sqrt{-h}N (Kf_\phi A_*-(\nabla_\nu n^\mu )n^\nu f_\phi\nabla_\mu \phi)$$

where I used identity $$K=\nabla_\nu n^\nu$$ and $$\nabla_\alpha f=f_\phi \nabla_\alpha \phi$$ (chain rule).

I have trouble with the last part of the above equation - i have no idea how to simplify this expression and put it into 3+1 form. Here are some references that I'm trying to follow: https://arxiv.org/abs/1101.3403, https://arxiv.org/abs/1708.02951, https://arxiv.org/abs/1812.02667, https://arxiv.org/abs/1512.06820.

• Minor comment to the post (v3): The word Langrangian should be Lagrangian. – Qmechanic Mar 7 at 18:17