Commit b9909ca9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

add a notebook just for fun

parent 3d4b1b18
Pipeline #23919 passed with stage
in 4 minutes and 42 seconds
%% Cell type:markdown id: tags:
# NIFTy4 tutorial
Import the necessary packages:
%% Cell type:code id: tags:
``` python
# some voodoo to make the notebook run properly
%matplotlib inline
import nifty4 as ift
import numpy as np
```
%% Cell type:markdown id: tags:
Define a space for our data
%% Cell type:code id: tags:
``` python
spc = ift.RGSpace(10) # 1D Cartesian space, contains 10 points
```
%% Cell type:markdown id: tags:
ask the space about its configuration:
%% Cell type:code id: tags:
``` python
print(spc)
```
%% Cell type:markdown id: tags:
Now, create a Field with Gaussian random numbers, living on that space
%% Cell type:code id: tags:
``` python
fld = ift.Field.from_random("normal",spc)
```
%% Cell type:code id: tags:
``` python
ift.plot(fld)
```
%% Cell type:markdown id: tags:
Now, let's do the same thing in two dimensions:
%% Cell type:code id: tags:
``` python
spc = ift.RGSpace((10,10))
fld = ift.Field.from_random("normal",spc)
ift.plot(fld)
```
%% Cell type:markdown id: tags:
Now, let's do some basic arithmetics with the field
%% Cell type:code id: tags:
``` python
fld += 50
fld *= -2
ift.plot(fld)
```
%% Cell type:markdown id: tags:
... and smooth it.
%% Cell type:code id: tags:
``` python
smooth = ift.FFTSmoothingOperator(spc,sigma=0.1)
fld_smooth = smooth(fld)
ift.plot(fld_smooth)
```
%% Cell type:code id: tags:
``` python
np.random.seed(43)
# Set up physical constants
# Total length of interval or volume the field lives on, e.g. in meters
L = 2.
# Typical distance over which the field is correlated (in same unit as L)
correlation_length = 0.1
# Variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance = 2.
# Smoothing length of response (in same unit as L)
response_sigma = 0.1
# Define resolution (pixels per dimension)
N_pixels = 256
# Set up derived constants
k_0 = 1./correlation_length
# Note that field_variance**2 = a*k_0/4. for this analytic form of power
# spectrum
a = field_variance**2/k_0*4.
pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
pixel_width = L/N_pixels
# Set up the geometry
s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width)
fft = ift.FFTOperator(s_space)
h_space = fft.target[0]
p_space = ift.PowerSpace(h_space)
# Create mock data
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
sp = ift.PS_field(p_space, pow_spec)
sh = ift.power_synthesize(sp, real_signal=True)
ss = fft.inverse_times(sh)
R = ift.FFTSmoothingOperator(s_space, sigma=response_sigma)
signal_to_noise = 1
diag = ift.Field(s_space, ss.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(diag)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
mean=0)
d = R(ss) + n
# Wiener filter
j = R.adjoint_times(N.inverse_times(d))
IC = ift.GradientNormController(verbose=True, iteration_limit=500,
tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC)
S_inv = fft.adjoint*Sh.inverse*fft
D = (R.adjoint*N.inverse*R + S_inv).inverse
# MR FIXME: we can/should provide a preconditioner here as well!
D = ift.InversionEnabler(D, inverter)
m = D(j)
```
%% Cell type:code id: tags:
``` python
ift.plot(m)
```
%% Cell type:code id: tags:
``` python
ift.plot(d)
```
%% Cell type:code id: tags:
``` python
ift.plot(ss)
```
%% Cell type:code id: tags:
``` python
```
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