Commit fc53aac7 authored by Vincent Eberle's avatar Vincent Eberle
Browse files

rng numpy update / nifty methods update / some small changes

parent 09de3f11
Pipeline #79930 passed with stage
in 36 seconds
......@@ -54,6 +54,7 @@
\maketitle
\VerbatimFootnotes
\section*{Nomenclature}
$d$ is the data vector and $s$ is the physical field about which one wants to
learn something given the data $d$ in a Bayesian fashion.
......@@ -67,11 +68,10 @@ building blocks.
the data, the statistics of the data including the error bars and a
description of the measurement device itself.
\item The \emph{prior} $\mathcal P(s)$ describes the knowledge the scientist has
\emph{before} executing the experiment. In order to define prior one needs to
\emph{before} executing the experiment. In order to define a prior one needs to
make one's knowledge about the physical process which was observed by the
measurement device explicit.
\item Finally, one needs an algorithmic and computational framework which is
able to actually do the Bayesian inference given the above information.
\item Finally, one needs an algorithmic and computational framework which is able to actually do the Bayesian inference, compute the \emph{posterior} $\mathcal P(s/d)$, given the above information.
\end{enumerate}
It becomes clear that these three parts are separate from each other. To
......@@ -83,8 +83,8 @@ The IFT group led by Torsten En{\ss}lin has knowledge and GPL-licensed software
to address the third part.
The interesting part is the second one, i.e. thinking about the priors. Here one
needs to have an understanding both of the underlying physics and probability
theory and Bayesian statistics. This is where an interesting interaction between
needs to have an understanding both of the underlying physics and Bayesian probability
theory. This is where an interesting interaction between
observers, i.e. you, and information field theorists, i.e. Torsten and his
group, will happen.
......@@ -121,10 +121,11 @@ That's it. The rest of this paper explains what these formulae mean and how to
compute them. From now one, our discussion will become increasingly specific.
\section*{Becoming more specific}
If the likelihood is Gaussian
If the likelihood is Gaussian (with noise covariance $N$ and response $R$)
\begin{align*}
\mathcal H(d|s) \propto (d-R(s))^\dagger N^{-1}(d-R(s))
\end{align*}
or Poissonian
\begin{align*}
\mathcal H(d|s) \propto - \log (R(s))^\dagger d+\sum_i R(s)_i,
......@@ -155,7 +156,7 @@ multi-frequency image of sky, etc. To cut a long story short: be aware of the
shape of \texttt{s} and \texttt{d}!
\subsection*{The response}
When we talk about \enquote{the response} we mean a function \texttt{R(s)} which
When we talk about \enquote{the response}, we mean a function \texttt{R(s)} which
takes an array of shape \texttt{s.shape} and returns an array of shape
\texttt{d.shape}. This function shall be made such that it simulates the
measurement device. Assume one knows the signal \texttt{s}, what would be the
......@@ -189,7 +190,7 @@ 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=\text{position}$ the gradient is a matrix with
\texttt{s.size} columns and \texttt{d.size} rows. Obviously, it is not
feasible to store such a matrix on a computer due to its size. There we think
feasible to store such a matrix on a computer due to its size. Therefore we think
of this matrix in terms of a linear map which maps an array of shape
\texttt{s.shape} to an array of shape \texttt{d.shape}. This linear map shall
be implemented on the computer in terms of a function. Think of this map as a
......@@ -214,17 +215,20 @@ we've got a linear function again which shall be the adjoint to
In order to test your implementation make sure that the following code runs
through:
\begin{lstlisting}
from numpy.random import default_rng
rng = default_rng()
# Test adjointness
position = np.random.random(s.shape)
s0 = np.random.random(s.shape)
d0 = np.random.random(d.shape)
position = rng.random(s.shape)
s0 = rng.random(s.shape)
d0 = rng.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 = rng.random(s.shape)
dir /= dir.var()
eps = 1e-8
finite_diff = (R(position + eps*dir)-R(position)) / eps
......@@ -240,9 +244,14 @@ NIFTy comes in. NIFTy has the notion of \enquote{spaces}. In a nutshell, a NIFTy
space comes with the following properties:
\begin{itemize}
\item Topology of the space: NIFTy supports regular grid spaces of arbitrary
dimensions, two 2-sphere pixelations and unstructured spaces (which do not
\item Topology of the space supported by NIFTY:
\begin{itemize}
\item regular grid spaces of arbitrary dimensions (RGSpace)
\item two 2-sphere pixelations, namely HEALPix (HPSpace) and Gauss-Legendre (GLSpace)
\item unstructured spaces (UnstructuredDomain) (which do not
have a specific topology and do not have pixel sizes).
\end{itemize}
\item Size of a pixel. This is called \enquote{distances} in NIFTy.
\end{itemize}
......@@ -263,20 +272,20 @@ Some tips exclusively for you:
\item Please have a look at the method \texttt{weight()} of \texttt{Field}. With
it, you can easily multiply the field values with or divide by the volume of
each pixel.
\item The methods \texttt{from\_global\_data()} and \texttt{to\_global\_data()}
\item The methods \texttt{from\_raw()} and \texttt{val()}
convert a Field to a numpy array and vice versa.
\item Be aware of the fact that NIFTy fields are immutable. As soon as you
pass a numpy array into a NIFTy field with \texttt{from\_global\_data()}, the
pass a numpy array into a NIFTy field with \texttt{from\_raw()}, the
\enquote{lock} flag will be set and the numpy array is immutable
afterwards. \texttt{to\_global\_data()} returns an object reference to this
afterwards. \texttt{val()} returns an object reference to this
numpy array which is locked and cannot be modified. If you want to modify
it, you may want to obtain an unlocked copy via
\texttt{to\_global\_data\_rw()}.
\texttt{val\_rw()}.
\end{itemize}
The point is: You need to fiddle around until the following test passes:
\begin{lstlisting}
import nifty5 as ift
import nifty6 as ift
# op is the linear operator you have defined
ift.extra.consistency_check(op)
\end{lstlisting}
......@@ -304,9 +313,10 @@ Why is this already sufficient?
\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
\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}
\end{document}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment