\documentclass{scrartcl}
\usepackage{csquotes}
\usepackage{hyperref}
\usepackage{fancyvrb}
\usepackage{listings}
\usepackage{color}
\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}
\lstset{
backgroundcolor=\color{white}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}; should come as last argument
basicstyle=\footnotesize, % the size of the fonts that are used for the code
breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
breaklines=true, % sets automatic line breaking
captionpos=b, % sets the caption-position to bottom
commentstyle=\color{mygreen}, % comment style
deletekeywords={...}, % if you want to delete keywords from the given language
escapeinside={\%*}{*)}, % if you want to add LaTeX within your code
extendedchars=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
frame=single, % adds a frame around the code
keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
keywordstyle=\color{blue}, % keyword style
language=Python, % the language of the code
morekeywords={*,...}, % if you want to add more keywords to the set
numbers=left, % where to put the line-numbers; possible values are (none, left, right)
numbersep=5pt, % how far the line-numbers are from the code
numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers
rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
showstringspaces=false, % underline spaces within strings only
showtabs=false, % show tabs within strings adding particular underscores
stepnumber=2, % the step between two line-numbers. If it's 1, each line will be numbered
stringstyle=\color{mymauve}, % string literal style
tabsize=2, % sets default tabsize to 2 spaces
title=\lstname % show the filename of files included with \lstinputlisting; also try caption instead of title
}
\usepackage{graphicx}
\usepackage{amsmath, amssymb}
\graphicspath{{img/}}
\usepackage[backend=biber, sorting=none]{biblatex}
\addbibresource{bib.bib}
\title{IFT Collaboration Guide}
\begin{document}
\VerbatimFootnotes
\section*{The need for a guide}
\section*{Mathematical description of what NIFTy needs}
For first order minimization:
\begin{itemize}
\item $\mathcal H (d|s) = - \log \mathcal P (d|s) + \text{constant in }s$
\item $\mathcal H' = \frac{\partial}{\partial s} \mathcal H (d|s)$
\end{itemize}
For second order minimization additionally:
\begin{itemize}
\item $\langle \mathcal H' \mathcal H'^\dagger \rangle_{\mathcal P (d|s)}$.
\end{itemize}
Note, that for Gaussian, Poissonian and Bernoulli likelihoods this term doesn't
need to be calculated and implemented because NIFTy does it automatically.
\section*{Becoming more specific}
If the likelihood is Gaussian ($\mathcal H(d|s) \propto (d-R(s))^\dagger N^{-1}
(d-R(s))$) or Poissonian ($\mathcal H(d|s) \propto - \log (R(s))^\dagger
d+\sum_i R(s)_i$ ), NIFTy needs:
\begin{itemize}
\item $R$.
\item $R'^\dagger := (\frac{\partial}{\partial s} R(s))^\dagger$.
\item $R' = \frac{\partial}{\partial s} R(s)$. (only for 2nd order minimization)
\end{itemize}
\section*{Even more specific}
Since NIFTy is implemented in python and is based on numpy let us be as specific
as possible and talk about numpy arrays. In the end, $s$ and $d$ will be a numpy
array defined in a python script.
The array \verb|d| is created by a script which reads in the actual data which
drops out of an instrument. This array will remain constant throughout the IFT
algorithm since the data is given and never will be changed.
The array \verb|s| is obviously not constant; only its shape won't change in the
course of the algorithm. NIFTy will store here the reconstruction of the
physical field which the user wants to infer. \verb|s| would be a
one-dimensional array for a time-series, two-dimensional for a multi-frequency
time-series or for a single-frequency image of the sky, three-dimensional for a
multi-frequency image of sky, etc. To cut a long story short: be aware of the
shape of \verb|s| and \verb|d|!
Now, the response appears. When we talk about \enquote{the response} we mean a
function \verb|R(s)| which takes an array of shape \verb|s.shape| and returns an
array of shape \verb|d.shape|. This function shall be made such that it
simulates the measurement device. Assume one knows the signal \verb|s|, what
would be the data \verb|d| in a noiseless measurement? Make sure that the
following code runs:
\begin{lstlisting}
import numpy as np
# Load the data and store it in a numpy array called 'd'.
# Store the shape of s in 'shp'.
# Store the function which implements the response in 'R'.
response_out = R(np.ones(shp))
if response_out.shape == d.shape:
print('Yay!')
else:
raise ValueError('Output of response has not the correct shape.')
\end{lstlisting}
Next, $R'$ and $R'^\dagger$ needs to be implemented. Since $R$ is a
function:
\begin{align*}
R: \mathbb R^{s.shape} \to \mathbb R^{d.shape},
\end{align*}
its gradient at the position \verb|s=position| is a linear map of the following shape:
\footnote{There are various ways to think about derivatives and gradients of
multi-dimensional functions. A different view on gradients would be that at a
given point $s=s_0$ the gradient is a matrix with \verb|s.size| columns and
\verb|d.size| rows. Obviously, it is not feasible to store such a matrix on a
computer due to its size. There we think of this matrix in terms of a linear
map which maps an array of shape \verb|s.shape| to an array of shape
\verb|d.shape|. This linear map shall be implemented on the computer in terms
of a function. Think of this map as a linear approximation to \verb|R| based
at \verb|s_0|.}
\begin{align*}
\left. \frac{dR}{ds}\right|_{s=position} = R': \mathbb R^{s.shape} \to \mathbb R^{d.shape}
\end{align*}
What needs to be implemented is a function \verb|R_prime(position, s0)| which
takes the arguments \verb|position| (which is an array of shape \verb|s.shape|
and determines the position at which we want to calculate the derivative) and
the array \verb|s0| which shall be the derivative taken of.
\verb|R_prime| is nonlinear in \verb|position| in general and linear in
\verb|s0|. The output of \verb|R_prime| is of shape \verb|d.shape|.
Finally, we need to implement the function \verb|R_prime_adjoint(postition,
d0)|. It implements the adjoint action of the derivative: $R'^\dagger$. It takes
a \verb|position| of shape \verb|s.shape| as well. After the position is set
we've got a linear function again which shall be the adjoint to
\verb|R_prime(position, _)|. Thus, \verb|d0| has the same shape as the data
\verb|d| and the output of \verb|R_prime| has the shape \verb|s.shape|.
In order to test your implementation make sure that the following code runs
through:
\begin{lstlisting}
# Test adjointness
position = np.random.random(s.shape)
s0 = np.random.random(s.shape)
d0 = np.random.random(d.shape)
res0 = d0.vdot(R_prime(position, s0))
res1 = s0.vdot(R_prime_adjoint(position, d0))
np.testing.assert_allclose(res0, res1)
# Test derivative
dir = np.random.random(s.shape)
dir /= dir.var()
eps = 1e-8
finite_diff = (R(position + eps*dir)-R(position)) / eps
derivative = R_prime(position, dir)
np.testing.assert_allclose(finite_diff, derivative)
\end{lstlisting}
\section*{The same in NIFTy language}
\section*{Example $\gamma$-ray imaging}
The information a $\gamma$-ray astronomer would provide to the algorithm (in the
simplest case):
\begin{itemize}
\item Data has Poissonian statistics.
\item Two functions: \verb|R(s)| and \verb|R_adjoint(d)| where $R$ applies a convolution with
a Gaussian kernel of given size to \verb|s| and picks out the positions to which
there exists a data point in \verb|d|. \verb|R_adjoint(d)| implements $R^\dagger$.
\end{itemize}
Why is this already sufficient?
\begin{itemize}
\item The Hamiltonian $\mathcal H$ is given by: $\mathcal H (d|s) = - \log
(R(s))^\dagger d+\sum_i R(s)_i$. Implementing $R$ and stating that the data is
Poissonian determines this form.
\item Since $R$ is a composition of a convolution and a sampling both of which
is a linear operation, $R$ itself is a linear operator.\footnote{I.e. $R(\alpha
s_1+s_2) = \alpha R(s_1) + R(s_2)$.} Thus, $R' = R$ and $R'^\dagger =
R^\dagger$. All in all, we need an implementation for $R$ and $R^\dagger$.
\end{itemize}
% \section*{Bibliography test}
% RESOLVE was first presented in \cite{Resolve2016}.
% \printbibliography
\end{document}