\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}