Commit f4d88ab7 authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'vincent_update' into 'master'

rng numpy update / nifty methods update / some small changes

See merge request !3
parents 269a8536 53d1c9b2
Pipeline #79998 passed with stage
in 30 seconds
......@@ -54,6 +54,7 @@
$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.
......@@ -70,8 +71,7 @@ building blocks.
\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.
It becomes clear that these three parts are separate from each other. To
......@@ -121,11 +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$)
\mathcal H(d|s) \propto (d-R(s))^\dagger N^{-1}(d-R(s))
with noise covariance $N$ and response $R$,\\ or Poissonian
or Poissonian
\mathcal H(d|s) \propto - \log (R(s))^\dagger d+\sum_i R(s)_i,
......@@ -214,17 +214,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
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 +243,14 @@ NIFTy comes in. NIFTy has the notion of \enquote{spaces}. In a nutshell, a NIFTy
space comes with the following properties:
\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:
\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).
\item Size of a pixel. This is called \enquote{distances} in NIFTy.
......@@ -263,20 +271,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
The point is: You need to fiddle around until the following test passes:
import nifty5 as ift
import nifty6 as ift
# op is the linear operator you have defined
......@@ -309,4 +317,5 @@ Why is this already sufficient?
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$.
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