Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
On Thursday, 7th July from 1 to 3 pm there will be a maintenance with a short downtime of GitLab.
Open sidebar
ift
NIFTy
Commits
1b7e9fb4
Commit
1b7e9fb4
authored
Nov 14, 2017
by
Martin Reinecke
Browse files
merge cleanups from dobj_experiments
parent
c7beff09
Pipeline
#21556
passed with stage
in 4 minutes and 50 seconds
Changes
10
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
demos/critical_filtering.py
View file @
1b7e9fb4
...
...
@@ -59,7 +59,7 @@ if __name__ == "__main__":
S
=
ift
.
create_power_operator
(
h_space
,
power_spectrum
=
p_spec
)
# Draw a sample sh from the prior distribution in harmonic space
sp
=
ift
.
F
ield
(
p_space
,
val
=
p_spec
(
p_space
.
k_lengths
)
)
sp
=
ift
.
PS_f
ield
(
p_space
,
p_spec
)
sh
=
ift
.
power_synthesize
(
sp
,
real_signal
=
True
)
# Choose the measurement instrument
...
...
@@ -106,7 +106,9 @@ if __name__ == "__main__":
def
ps0
(
k
):
return
(
1.
/
(
1.
+
k
)
**
2
)
t0
=
ift
.
Field
(
p_space
,
val
=
np
.
log
(
1.
/
(
1
+
p_space
.
k_lengths
)
**
2
))
t0
=
ift
.
Field
(
p_space
,
val
=
ift
.
dobj
.
from_global_data
(
np
.
log
(
1.
/
(
1
+
p_space
.
k_lengths
)
**
2
)))
for
i
in
range
(
500
):
S0
=
ift
.
create_power_operator
(
h_space
,
power_spectrum
=
ps0
)
...
...
@@ -138,6 +140,6 @@ if __name__ == "__main__":
t0
=
power_energy
.
position
.
real
# Plot current estimate
print
(
i
)
ift
.
dobj
.
m
print
(
i
)
if
i
%
5
==
0
:
plot_parameters
(
m0
,
t0
,
ift
.
log
(
sp
),
data_power
)
demos/log_normal_wiener_filter.py
View file @
1b7e9fb4
...
...
@@ -3,35 +3,38 @@ import numpy as np
if
__name__
==
"__main__"
:
np
.
random
.
seed
(
42
)
# Setting up parameters
|\label{code:wf_parameters}|
# Setting up parameters
correlation_length
=
1.
# Typical distance over which the field is correlated
field_variance
=
2.
# Variance of field in position space
response_sigma
=
0.02
# Smoothing length of response (in same unit as L)
signal_to_noise
=
100
# The signal to noise ratio
np
.
random
.
seed
(
43
)
# Fixing the random seed
def
power_spectrum
(
k
):
# Defining the power spectrum
a
=
4
*
correlation_length
*
field_variance
**
2
return
a
/
(
1
+
k
*
correlation_length
)
**
4
# Setting up the geometry
|\label{code:wf_geometry}|
L
=
2.
# Total side-length of the domain
N_pixels
=
128
# Grid resolution (pixels per axis)
#signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
# Setting up the geometry
L
=
2.
# Total side-length of the domain
N_pixels
=
128
# Grid resolution (pixels per axis)
#
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space
=
ift
.
HPSpace
(
16
)
harmonic_space
=
signal_space
.
get_default_codomain
()
fft
=
ift
.
FFTOperator
(
harmonic_space
,
target
=
signal_space
)
power_space
=
ift
.
PowerSpace
(
harmonic_space
)
# Creating the mock signal |\label{code:wf_mock_signal}|
S
=
ift
.
create_power_operator
(
harmonic_space
,
power_spectrum
=
power_spectrum
)
mock_power
=
ift
.
Field
(
power_space
,
val
=
power_spectrum
(
power_space
.
k_lengths
))
# Creating the mock signal
S
=
ift
.
create_power_operator
(
harmonic_space
,
power_spectrum
=
power_spectrum
)
mock_power
=
ift
.
PS_field
(
power_space
,
power_spectrum
)
mock_signal
=
fft
(
ift
.
power_synthesize
(
mock_power
,
real_signal
=
True
))
# Setting up an exemplary response
mask
=
ift
.
Field
.
ones
(
signal_space
)
N10
=
int
(
N_pixels
/
10
)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R
=
ift
.
ResponseOperator
(
signal_space
,
sigma
=
(
response_sigma
,),
exposure
=
(
mask
,))
#|\label{code:wf_response}|
# mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R
=
ift
.
ResponseOperator
(
signal_space
,
sigma
=
(
response_sigma
,),
exposure
=
(
mask
,))
data_domain
=
R
.
target
[
0
]
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
])
...
...
@@ -39,40 +42,52 @@ if __name__ == "__main__":
ndiag
=
ift
.
Field
.
full
(
data_domain
,
mock_signal
.
var
()
/
signal_to_noise
)
N
=
ift
.
DiagonalOperator
(
ndiag
.
weight
(
1
))
noise
=
ift
.
Field
.
from_random
(
domain
=
data_domain
,
random_type
=
'normal'
,
std
=
mock_signal
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
data
=
R
(
ift
.
exp
(
mock_signal
))
+
noise
#|\label{code:wf_mock_data}|
std
=
mock_signal
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
data
=
R
(
ift
.
exp
(
mock_signal
))
+
noise
# Wiener filter
m0
=
ift
.
Field
.
zeros
(
harmonic_space
)
ctrl
=
ift
.
GradientNormController
(
verbose
=
False
,
tol_abs_gradnorm
=
1
)
ctrl2
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
0.1
,
name
=
"outer"
)
ctrl
=
ift
.
GradientNormController
(
verbose
=
False
,
tol_abs_gradnorm
=
1
)
ctrl2
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
0.1
,
name
=
"outer"
)
inverter
=
ift
.
ConjugateGradient
(
controller
=
ctrl
)
energy
=
ift
.
library
.
LogNormalWienerFilterEnergy
(
m0
,
data
,
R_harmonic
,
N
,
S
,
inverter
=
inverter
)
minimizer1
=
ift
.
VL_BFGS
(
controller
=
ctrl2
,
max_history_length
=
20
)
energy
=
ift
.
library
.
LogNormalWienerFilterEnergy
(
m0
,
data
,
R_harmonic
,
N
,
S
,
inverter
=
inverter
)
# minimizer1 = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer2
=
ift
.
RelaxedNewton
(
controller
=
ctrl2
)
minimizer3
=
ift
.
SteepestDescent
(
controller
=
ctrl2
)
#
minimizer3 = ift.SteepestDescent(controller=ctrl2)
#me1 = minimizer1(energy)
#
me1 = minimizer1(energy)
me2
=
minimizer2
(
energy
)
#me3 = minimizer3(energy)
#
me3 = minimizer3(energy)
#m1 = fft(me1[0].position)
#
m1 = fft(me1[0].position)
m2
=
fft
(
me2
[
0
].
position
)
#m3 = fft(me3[0].position)
# m3 = fft(me3[0].position)
#Plotting #|\label{code:wf_plotting}|
ift
.
plotting
.
plot
(
mock_signal
.
real
,
name
=
'mock_signal.pdf'
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
ift
.
plotting
.
plot
(
ift
.
Field
(
signal_space
,
val
=
np
.
log
(
data
.
val
.
real
).
reshape
(
signal_space
.
shape
)),
name
=
"log_of_data.pdf"
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
#ift.plotting.plot(m1.real,name='m_LBFGS.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index")
ift
.
plotting
.
plot
(
m2
.
real
,
name
=
'm_Newton.pdf'
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
#ift.plotting.plot(m3.real,name='m_SteepestDescent.pdf', colormap="plasma",xlabel="Pixel Index",ylabel="Pixel Index")
# Plotting
ift
.
plotting
.
plot
(
mock_signal
.
real
,
name
=
'mock_signal.pdf'
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
logdata
=
np
.
log
(
ift
.
dobj
.
to_global_data
(
data
.
val
.
real
)).
reshape
(
signal_space
.
shape
)
ift
.
plotting
.
plot
(
ift
.
Field
(
signal_space
,
val
=
ift
.
dobj
.
from_global_data
(
logdata
)),
name
=
"log_of_data.pdf"
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
# ift.plotting.plot(m1.real,name='m_LBFGS.pdf', colormap="plasma",
# xlabel="Pixel Index", ylabel="Pixel Index")
ift
.
plotting
.
plot
(
m2
.
real
,
name
=
'm_Newton.pdf'
,
colormap
=
"plasma"
,
xlabel
=
"Pixel Index"
,
ylabel
=
"Pixel Index"
)
# ift.plotting.plot(m3.real, name='m_SteepestDescent.pdf',
# colormap="plasma", xlabel="Pixel Index",
# ylabel="Pixel Index")
# Probing the variance
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
proby
=
Proby
(
signal_space
,
probe_count
=
1
)
proby
(
lambda
z
:
fft
(
me2
[
0
].
curvature
.
inverse_times
(
fft
.
adjoint_times
(
z
))))
sm
=
ift
.
FFTSmoothingOperator
(
signal_space
,
sigma
=
0.02
)
variance
=
sm
(
proby
.
diagonal
.
weight
(
-
1
))
ift
.
plotting
.
plot
(
variance
,
name
=
'variance.pdf'
)
ift
.
plotting
.
plot
(
variance
,
name
=
'variance.pdf'
)
demos/paper_demos/cartesian_wiener_filter.py
View file @
1b7e9fb4
...
...
@@ -2,82 +2,82 @@ import numpy as np
import
nifty2go
as
ift
if
__name__
==
"__main__"
:
signal_to_noise
=
1.5
# The signal to noise ratio
signal_to_noise
=
1.5
# The signal to noise ratio
# Setting up parameters
|\label{code:wf_parameters}|
correlation_length_1
=
1.
# Typical distance over which the field is correlated
field_variance_1
=
2.
# Variance of field in position space
# Setting up parameters
correlation_length_1
=
1.
# Typical distance over which the field is correlated
field_variance_1
=
2.
# Variance of field in position space
response_sigma_1
=
0.05
# Smoothing length of response (in same unit as L)
response_sigma_1
=
0.05
# Smoothing length of response (in same unit as L)
def
power_spectrum_1
(
k
):
# note: field_variance**2 = a*k_0/4.
def
power_spectrum_1
(
k
):
# note: field_variance**2 = a*k_0/4.
a
=
4
*
correlation_length_1
*
field_variance_1
**
2
return
a
/
(
1
+
k
*
correlation_length_1
)
**
4.
# Setting up the geometry
|\label{code:wf_geometry}|
L_1
=
2.
# Total side-length of the domain
N_pixels_1
=
512
# Grid resolution (pixels per axis)
# Setting up the geometry
L_1
=
2.
# Total side-length of the domain
N_pixels_1
=
512
# Grid resolution (pixels per axis)
signal_space_1
=
ift
.
RGSpace
([
N_pixels_1
],
distances
=
L_1
/
N_pixels_1
)
harmonic_space_1
=
signal_space_1
.
get_default_codomain
()
# Setting up the geometry
|\label{code:wf_geometry}|
L_2
=
2.
# Total side-length of the domain
N_pixels_2
=
512
# Grid resolution (pixels per axis)
# Setting up the geometry
L_2
=
2.
# Total side-length of the domain
N_pixels_2
=
512
# Grid resolution (pixels per axis)
signal_space_2
=
ift
.
RGSpace
([
N_pixels_2
],
distances
=
L_2
/
N_pixels_2
)
harmonic_space_2
=
signal_space_2
.
get_default_codomain
()
signal_domain
=
ift
.
DomainTuple
.
make
((
signal_space_1
,
signal_space_2
))
mid_domain
=
ift
.
DomainTuple
.
make
((
signal_space_1
,
harmonic_space_2
))
harmonic_domain
=
ift
.
DomainTuple
.
make
((
harmonic_space_1
,
harmonic_space_2
))
harmonic_domain
=
ift
.
DomainTuple
.
make
((
harmonic_space_1
,
harmonic_space_2
))
fft_1
=
ift
.
FFTOperator
(
harmonic_domain
,
space
=
0
)
power_space_1
=
ift
.
PowerSpace
(
harmonic_space_1
)
mock_power_1
=
ift
.
F
ield
(
power_space_1
,
val
=
power_spectrum_1
(
power_space_1
.
k_lengths
)
)
mock_power_1
=
ift
.
PS_f
ield
(
power_space_1
,
power_spectrum_1
)
# Setting up parameters
correlation_length_2
=
1.
# Typical distance over which the field is correlated
field_variance_2
=
2.
# Variance of field in position space
response_sigma_2
=
0.01
# Smoothing length of response (in same unit as L)
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_2
=
1.
# Typical distance over which the field is correlated
field_variance_2
=
2.
# Variance of field in position space
response_sigma_2
=
0.01
# Smoothing length of response (in same unit as L)
def
power_spectrum_2
(
k
):
# note: field_variance**2 = a*k_0/4.
def
power_spectrum_2
(
k
):
# note: field_variance**2 = a*k_0/4.
a
=
4
*
correlation_length_2
*
field_variance_2
**
2
return
a
/
(
1
+
k
*
correlation_length_2
)
**
2.5
fft_2
=
ift
.
FFTOperator
(
mid_domain
,
space
=
1
)
power_space_2
=
ift
.
PowerSpace
(
harmonic_space_2
)
mock_power_2
=
ift
.
F
ield
(
power_space_2
,
val
=
power_spectrum_2
(
power_space_2
.
k_lengths
)
)
mock_power_2
=
ift
.
PS_f
ield
(
power_space_2
,
power_spectrum_2
)
fft
=
ift
.
ComposedOperator
((
fft_1
,
fft_2
))
mock_power
=
ift
.
Field
(
domain
=
(
power_space_1
,
power_space_2
),
val
=
np
.
outer
(
mock_power_1
.
val
,
mock_power_2
.
val
))
val
=
ift
.
dobj
.
from_global_data
(
np
.
outer
(
ift
.
dobj
.
to_global_data
(
mock_power_1
.
val
)
,
ift
.
dobj
.
to_global_data
(
mock_power_2
.
val
))
))
diagonal
=
ift
.
power_synthesize_special
(
mock_power
,
spaces
=
(
0
,
1
))
**
2
diagonal
=
diagonal
.
real
S
=
ift
.
DiagonalOperator
(
diagonal
)
np
.
random
.
seed
(
10
)
mock_signal
=
fft
(
ift
.
power_synthesize
(
mock_power
,
real_signal
=
True
))
# Setting up a exemplary response
N1_10
=
int
(
N_pixels_1
/
10
)
mask_1
=
ift
.
Field
.
ones
(
signal_space_1
)
mask_1
.
val
[
N1_10
*
7
:
N1_10
*
9
]
=
0.
mask_1
=
np
.
ones
(
signal_space_1
.
shape
)
mask_1
[
N1_10
*
7
:
N1_10
*
9
]
=
0.
mask_1
=
ift
.
Field
(
signal_space_1
,
ift
.
dobj
.
from_global_data
(
mask_1
))
N2_10
=
int
(
N_pixels_2
/
10
)
mask_2
=
ift
.
Field
.
ones
(
signal_space_2
)
mask_2
.
val
[
N2_10
*
7
:
N2_10
*
9
]
=
0.
mask_2
=
np
.
ones
(
signal_space_2
.
shape
)
mask_2
[
N2_10
*
7
:
N2_10
*
9
]
=
0.
mask_2
=
ift
.
Field
(
signal_space_2
,
ift
.
dobj
.
from_global_data
(
mask_2
))
R
=
ift
.
ResponseOperator
(
signal_domain
,
spaces
=
(
0
,
1
),
R
=
ift
.
ResponseOperator
(
signal_domain
,
spaces
=
(
0
,
1
),
sigma
=
(
response_sigma_1
,
response_sigma_2
),
exposure
=
(
mask_1
,
mask_2
))
#|\label{code:wf_response}|
exposure
=
(
mask_1
,
mask_2
))
data_domain
=
R
.
target
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
])
...
...
@@ -87,21 +87,23 @@ if __name__ == "__main__":
noise
=
ift
.
Field
.
from_random
(
domain
=
data_domain
,
random_type
=
'normal'
,
std
=
mock_signal
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
data
=
R
(
mock_signal
)
+
noise
#|\label{code:wf_mock_data}|
data
=
R
(
mock_signal
)
+
noise
# Wiener filter
j
=
R_harmonic
.
adjoint_times
(
N
.
inverse_times
(
data
))
ctrl
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
0.1
)
inverter
=
ift
.
ConjugateGradient
(
controller
=
ctrl
)
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
wiener_curvature
=
ift
.
InversionEnabler
(
wiener_curvature
,
inverter
)
m_k
=
wiener_curvature
.
inverse_times
(
j
)
#|\label{code:wf_wiener_filter}|
m_k
=
wiener_curvature
.
inverse_times
(
j
)
m
=
fft
(
m_k
)
# Probing the variance
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
proby
=
Proby
((
signal_space_1
,
signal_space_2
),
probe_count
=
1
,
ncpu
=
1
)
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
proby
=
Proby
((
signal_space_1
,
signal_space_2
),
probe_count
=
1
,
ncpu
=
1
)
proby
(
lambda
z
:
fft
(
wiener_curvature
.
inverse_times
(
fft
.
inverse_times
(
z
))))
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
...
...
demos/paper_demos/wiener_filter.py
View file @
1b7e9fb4
...
...
@@ -3,34 +3,38 @@ import numpy as np
if
__name__
==
"__main__"
:
# Setting up parameters
|\label{code:wf_parameters}|
# Setting up parameters
correlation_length_scale
=
1.
# Typical distance over which the field is correlated
fluctuation_scale
=
2.
# Variance of field in position space
response_sigma
=
0.05
# Smoothing length of response (in same unit as L)
signal_to_noise
=
1.5
# The signal to noise ratio
np
.
random
.
seed
(
43
)
# Fixing the random seed
def
power_spectrum
(
k
):
# Defining the power spectrum
a
=
4
*
correlation_length_scale
*
fluctuation_scale
**
2
return
a
/
(
1
+
(
k
*
correlation_length_scale
)
**
2
)
**
2
# Setting up the geometry
|\label{code:wf_geometry}|
# Setting up the geometry
L
=
2.
# Total side-length of the domain
N_pixels
=
512
# Grid resolution (pixels per axis)
signal_space
=
ift
.
RGSpace
([
N_pixels
,
N_pixels
],
distances
=
L
/
N_pixels
)
harmonic_space
=
signal_space
.
get_default_codomain
()
fft
=
ift
.
FFTOperator
(
harmonic_space
,
target
=
signal_space
)
power_space
=
ift
.
PowerSpace
(
harmonic_space
,
binbounds
=
ift
.
PowerSpace
.
useful_binbounds
(
harmonic_space
,
logarithmic
=
True
))
power_space
=
ift
.
PowerSpace
(
harmonic_space
,
binbounds
=
ift
.
PowerSpace
.
useful_binbounds
(
harmonic_space
,
logarithmic
=
True
))
# Creating the mock signal |\label{code:wf_mock_signal}|
S
=
ift
.
create_power_operator
(
harmonic_space
,
power_spectrum
=
power_spectrum
)
mock_power
=
ift
.
Field
(
power_space
,
val
=
power_spectrum
(
power_space
.
k_lengths
))
# Creating the mock signal
S
=
ift
.
create_power_operator
(
harmonic_space
,
power_spectrum
=
power_spectrum
)
mock_power
=
ift
.
PS_field
(
power_space
,
power_spectrum
)
mock_signal
=
fft
(
ift
.
power_synthesize
(
mock_power
,
real_signal
=
True
))
# Setting up an exemplary response
mask
=
ift
.
Field
.
ones
(
signal_space
)
mask
=
np
.
ones
(
signal_space
.
shape
)
N10
=
int
(
N_pixels
/
10
)
mask
.
val
[
N10
*
5
:
N10
*
9
,
N10
*
5
:
N10
*
9
]
=
0.
R
=
ift
.
ResponseOperator
(
signal_space
,
sigma
=
(
response_sigma
,),
exposure
=
(
mask
,))
#|\label{code:wf_response}|
mask
[
N10
*
5
:
N10
*
9
,
N10
*
5
:
N10
*
9
]
=
0.
mask
=
ift
.
Field
(
signal_space
,
ift
.
dobj
.
from_global_data
(
mask
))
R
=
ift
.
ResponseOperator
(
signal_space
,
sigma
=
(
response_sigma
,),
exposure
=
(
mask
,))
data_domain
=
R
.
target
[
0
]
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
])
...
...
@@ -39,28 +43,29 @@ if __name__ == "__main__":
N
=
ift
.
DiagonalOperator
(
ndiag
.
weight
(
1
))
noise
=
ift
.
Field
.
from_random
(
domain
=
data_domain
,
random_type
=
'normal'
,
std
=
mock_signal
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
data
=
R
(
mock_signal
)
+
noise
#|\label{code:wf_mock_data}|
data
=
R
(
mock_signal
)
+
noise
# Wiener filter
j
=
R_harmonic
.
adjoint_times
(
N
.
inverse_times
(
data
))
ctrl
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
0.1
)
ctrl
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
0.1
)
inverter
=
ift
.
ConjugateGradient
(
controller
=
ctrl
)
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
wiener_curvature
=
ift
.
InversionEnabler
(
wiener_curvature
,
inverter
)
m_k
=
wiener_curvature
.
inverse_times
(
j
)
#|\label{code:wf_wiener_filter}|
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
wiener_curvature
=
ift
.
InversionEnabler
(
wiener_curvature
,
inverter
)
m_k
=
wiener_curvature
.
inverse_times
(
j
)
m
=
fft
(
m_k
)
# Probing the uncertainty |\label{code:wf_uncertainty_probing}|
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
proby
=
Proby
(
signal_space
,
probe_count
=
1
,
ncpu
=
1
)
proby
(
lambda
z
:
fft
(
wiener_curvature
.
inverse_times
(
fft
.
inverse_times
(
z
))))
#|\label{code:wf_variance_fft_wrap}|
# Probing the uncertainty
class
Proby
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
proby
=
Proby
(
signal_space
,
probe_count
=
1
,
ncpu
=
1
)
proby
(
lambda
z
:
fft
(
wiener_curvature
.
inverse_times
(
fft
.
inverse_times
(
z
))))
sm
=
ift
.
FFTSmoothingOperator
(
signal_space
,
sigma
=
0.03
)
variance
=
ift
.
sqrt
(
sm
(
proby
.
diagonal
.
weight
(
-
1
)))
#|\label{code:wf_variance_weighting}|
variance
=
ift
.
sqrt
(
sm
(
proby
.
diagonal
.
weight
(
-
1
)))
# Plotting
#|\label{code:wf_plotting}|
# Plotting
ift
.
plotting
.
plot
(
variance
,
name
=
"uncertainty.pdf"
,
xlabel
=
'Pixel index'
,
ylabel
=
'Pixel index'
)
ift
.
plotting
.
plot
(
mock_signal
,
name
=
"mock_signal.pdf"
,
xlabel
=
'Pixel index'
,
ylabel
=
'Pixel index'
)
ift
.
plotting
.
plot
(
ift
.
Field
(
signal_space
,
val
=
data
.
val
),
name
=
"data.pdf"
,
xlabel
=
'Pixel index'
,
ylabel
=
'Pixel index'
)
ift
.
plotting
.
plot
(
m
,
name
=
"map.pdf"
,
xlabel
=
'Pixel index'
,
ylabel
=
'Pixel index'
)
demos/probing.py
View file @
1b7e9fb4
from
__future__
import
print_function
import
nifty2go
as
ift
import
numpy
as
np
np
.
random
.
seed
(
42
)
class
DiagonalProber
(
ift
.
DiagonalProberMixin
,
ift
.
Prober
):
pass
...
...
@@ -19,9 +19,9 @@ diagOp = ift.DiagonalOperator(f)
diagProber
=
DiagonalProber
(
domain
=
x
)
diagProber
(
diagOp
)
print
((
f
-
diagProber
.
diagonal
).
norm
())
ift
.
dobj
.
m
print
((
f
-
diagProber
.
diagonal
).
norm
())
multiProber
=
MultiProber
(
domain
=
x
)
multiProber
(
diagOp
)
print
((
f
-
multiProber
.
diagonal
).
norm
())
print
(
f
.
sum
()
-
multiProber
.
trace
)
ift
.
dobj
.
m
print
((
f
-
multiProber
.
diagonal
).
norm
())
ift
.
dobj
.
m
print
(
f
.
sum
()
-
multiProber
.
trace
)
demos/wiener_filter_via_curvature.py
View file @
1b7e9fb4
...
...
@@ -3,6 +3,8 @@ import nifty2go as ift
import
numericalunits
as
nu
if
__name__
==
"__main__"
:
# In MPI mode, the random seed for numericalunits must be set by hand
nu
.
reset_units
(
43
)
dimensionality
=
2
np
.
random
.
seed
(
43
)
...
...
@@ -43,8 +45,7 @@ if __name__ == "__main__":
power_spectrum
=
power_spectrum
)
np
.
random
.
seed
(
43
)
mock_power
=
ift
.
Field
(
power_space
,
val
=
power_spectrum
(
power_space
.
k_lengths
))
mock_power
=
ift
.
PS_field
(
power_space
,
power_spectrum
)
mock_harmonic
=
ift
.
power_synthesize
(
mock_power
,
real_signal
=
True
)
mock_harmonic
=
mock_harmonic
.
real
mock_signal
=
fft
(
mock_harmonic
)
...
...
@@ -80,7 +81,7 @@ if __name__ == "__main__":
ift
.
plotting
.
plot
(
ift
.
Field
(
sspace2
,
mock_signal
.
real
.
val
)
/
nu
.
K
,
name
=
"mock_signal.pdf"
)
ift
.
plotting
.
plot
(
ift
.
Field
(
sspace2
,
val
=
data
.
val
.
real
.
reshape
(
signal_space
.
shape
))
/
nu
.
K
,
name
=
"data.pdf"
)
data
=
ift
.
dobj
.
to_global_data
(
data
.
val
.
real
).
reshape
(
sspace2
.
shape
)
/
nu
.
K
data
=
ift
.
Field
(
sspace2
,
val
=
ift
.
dobj
.
from_global_data
(
data
))
/
nu
.
K
ift
.
plotting
.
plot
(
ift
.
Field
(
sspace2
,
val
=
data
),
name
=
"data.pdf"
)
ift
.
plotting
.
plot
(
ift
.
Field
(
sspace2
,
m_s
.
real
.
val
)
/
nu
.
K
,
name
=
"map.pdf"
)
demos/wiener_filter_via_hamiltonian.py
View file @
1b7e9fb4
...
...
@@ -50,14 +50,15 @@ if __name__ == "__main__":
S
=
ift
.
create_power_operator
(
h_space
,
power_spectrum
=
p_spec
)
# Drawing a sample sh from the prior distribution in harmonic space
sp
=
ift
.
F
ield
(
p_space
,
val
=
p_spec
(
p_space
.
k_lengths
)
)
sp
=
ift
.
PS_f
ield
(
p_space
,
p_spec
)
sh
=
ift
.
power_synthesize
(
sp
,
real_signal
=
True
)
ss
=
fft
.
adjoint_times
(
sh
)
# Choosing the measurement instrument
# Instrument = ift.FFTSmoothingOperator(s_space, sigma=0.05)
diag
=
ift
.
Field
.
ones
(
s_space
)
diag
.
val
[
20
:
80
,
20
:
80
]
=
0
diag
=
np
.
ones
(
s_space
.
shape
)
diag
[
20
:
80
,
20
:
80
]
=
0
diag
=
ift
.
Field
(
s_space
,
ift
.
dobj
.
from_global_data
(
diag
))
Instrument
=
ift
.
DiagonalOperator
(
diag
)
# Adding a harmonic transformation to the instrument
...
...
nifty/data_objects/numpy_do.py
View file @
1b7e9fb4
# Data object module for NIFTy that uses simple numpy ndarrays.
from
__future__
import
print_function
import
numpy
as
np
from
numpy
import
ndarray
as
data_object
from
numpy
import
full
,
empty
,
empty_like
,
sqrt
,
ones
,
zeros
,
vdot
,
abs
,
\
bincount
,
exp
,
log
exp
,
log
from
.random
import
Random
__all__
=
[
"ntask"
,
"rank"
,
"master"
,
"local_shape"
,
"data_object"
,
"full"
,
"empty"
,
"zeros"
,
"ones"
,
"empty_like"
,
"vdot"
,
"abs"
,
"exp"
,
"log"
,
"sqrt"
,
"bincount"
,
"from_object"
,
"from_random"
,
"local_data"
,
"ibegin"
,
"np_allreduce_sum"
,
"distaxis"
,
"from_local_data"
,
"from_global_data"
,
"to_global_data"
,
"redistribute"
,
"default_distaxis"
]
ntask
=
1
rank
=
0
master
=
True
def
mprint
(
*
args
):
print
(
*
args
)
def
from_object
(
object
,
dtype
=
None
,
copy
=
True
):
return
np
.
array
(
object
,
dtype
=
dtype
,
copy
=
copy
)
...
...
nifty/dobj.py
View file @
1b7e9fb4
from
.data_objects.numpy_do
import
*
__all__
=
[
"ntask"
,
"rank"
,
"master"
,
"local_shape"
,
"data_object"
,
"full"
,
"empty"
,
"zeros"
,
"ones"
,
"empty_like"
,
"vdot"
,
"abs"
,
"exp"
,
"log"
,
"sqrt"
,
"from_object"
,
"from_random"
,
"local_data"
,
"ibegin"
,
"np_allreduce_sum"
,
"distaxis"
,
"from_local_data"
,
"from_global_data"
,
"to_global_data"
,
"redistribute"
,
"default_distaxis"
,
"mprint"
]
nifty/sugar.py
View file @
1b7e9fb4
...
...
@@ -17,18 +17,12 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
import
numpy
as
np
from
.
import
Space
,
\
PowerSpace
,
\
Field
,
\
ComposedOperator
,
\
DiagonalOperator
,
\
PowerProjectionOperator
,
\
FFTOperator
,
\
sqrt
,
\
DomainTuple
from
.
import
Space
,
PowerSpace
,
Field
,
ComposedOperator
,
DiagonalOperator
,
\
PowerProjectionOperator
,
FFTOperator
,
sqrt
,
DomainTuple
,
dobj
from
.
import
nifty_utilities
as
utilities
__all__
=
[
'power_analyze'
,
__all__
=
[
'PS_field'
,
'power_analyze'
,
'power_synthesize'
,
'power_synthesize_special'
,
'create_power_field'
,
...
...
@@ -37,6 +31,13 @@ __all__ = ['power_analyze',
'create_composed_fft_operator'
]
def
PS_field
(
pspace
,
func
,
dtype
=
None
):
if
not
isinstance
(
pspace
,
PowerSpace
):
raise
TypeError
data
=
dobj
.
from_global_data
(
func
(
pspace
.
k_lengths
))
return
Field
(
pspace
,
val
=
data
,
dtype
=
dtype
)
def
_single_power_analyze
(
field
,
idx
,
binbounds
):
from
.operators.power_projection_operator
import
PowerProjectionOperator
power_domain
=
PowerSpace
(
field
.
domain
[
idx
],
binbounds
)
...
...
@@ -91,8 +92,8 @@ def power_analyze(field, spaces=None, binbounds=None,
# power_space instances
for
sp
in
field
.
domain
:
if
not
sp
.
harmonic
and
not
isinstance
(
sp
,
PowerSpace
):
print
(
"WARNING: Field has a space in `domain` which is "
"neither harmonic nor a PowerSpace."
)
dobj
.
m
print
(
"WARNING: Field has a space in `domain` which is "
"neither harmonic nor a PowerSpace."
)
# check if the `spaces` input is valid
if
spaces
is
None
:
...
...
@@ -222,8 +223,7 @@ def create_power_field(domain, power_spectrum, dtype=None):
fp
=
Field
(
power_domain
,
val
=
power_spectrum
.
val
,
dtype
=
dtype
)
else
:
power_domain
=
PowerSpace
(
domain
)
fp
=
Field
(
power_domain
,
val
=
power_spectrum
(
power_domain
.
k_lengths
),
dtype
=
dtype
)
fp
=
PS_field
(
power_domain
,
power_spectrum
,
dtype
)
P
=
PowerProjectionOperator
(
domain
,
power_domain
)
f
=
P
.
adjoint_times
(
fp
)
...
...
@@ -258,7 +258,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
"""
domain
=
DomainTuple
.
make
(
domain
)
if
space
is
None
:
if
len
(
domain
)
!=
1
:
if
len
(
domain
)
!=
1
:
raise
ValueError
(
"space keyword must be set"
)
else
:
space
=
0
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment