main.tex 8.94 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
4
\documentclass{scrartcl}

\usepackage{csquotes}
\usepackage{hyperref}
Philipp Arras's avatar
Philipp Arras committed
5
\usepackage{fancyvrb}
Philipp Arras's avatar
Philipp Arras committed
6

Philipp Arras's avatar
Philipp Arras committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
\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}
Philipp Arras's avatar
Philipp Arras committed
45
46
47
48
49
50
51
52
\graphicspath{{img/}}

\usepackage[backend=biber, sorting=none]{biblatex}
\addbibresource{bib.bib}

\title{IFT Collaboration Guide}

\begin{document}
Philipp Arras's avatar
Philipp Arras committed
53
\VerbatimFootnotes
Philipp Arras's avatar
Philipp Arras committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
\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}


Philipp Arras's avatar
Philipp Arras committed
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
\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}


Philipp Arras's avatar
Philipp Arras committed
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
\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}
Philipp Arras's avatar
Philipp Arras committed
195
196
197
198
199
200

% \section*{Bibliography test}
% RESOLVE was first presented in \cite{Resolve2016}.
% \printbibliography

\end{document}