main.tex 8.94 KB
 Philipp Arras committed Aug 24, 2018 1 2 3 4 \documentclass{scrartcl} \usepackage{csquotes} \usepackage{hyperref}  Philipp Arras committed Aug 24, 2018 5 \usepackage{fancyvrb}  Philipp Arras committed Aug 24, 2018 6   Philipp Arras committed Aug 24, 2018 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 committed Aug 24, 2018 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 committed Aug 24, 2018 53 \VerbatimFootnotes  Philipp Arras committed Aug 24, 2018 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 committed Aug 24, 2018 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 committed Aug 24, 2018 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 committed Aug 24, 2018 195 196 197 198 199 200  % \section*{Bibliography test} % RESOLVE was first presented in \cite{Resolve2016}. % \printbibliography \end{document}