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
Open sidebar
ift
NIFTy
Commits
07bcecbd
Commit
07bcecbd
authored
Nov 28, 2017
by
Martin Reinecke
Browse files
revert accidental commit
parent
8122ea7a
Pipeline
#22293
failed with stage
in 4 minutes and 1 second
Changes
1
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
demos/wiener_filter_via_curvature.py
View file @
07bcecbd
use_nifty2go
=
True
import
numpy
as
np
if
use_nifty2go
:
import
nifty2go
as
ift
else
:
import
nifty
as
ift
import
nifty2go
as
ift
import
numericalunits
as
nu
if
__name__
==
"__main__"
:
# In MPI mode, the random seed for numericalunits must be set by hand
if
not
use_nifty2go
:
ift
.
nifty_configuration
[
'default_distribution_strategy'
]
=
'fftw'
ift
.
nifty_configuration
[
'harmonic_rg_base'
]
=
'real'
nu
.
reset_units
(
43
)
dimensionality
=
2
np
.
random
.
seed
(
43
)
...
...
@@ -40,14 +32,11 @@ if __name__ == "__main__":
# Total side-length of the domain
L
=
2.
*
nu
.
m
# Grid resolution (pixels per axis)
N_pixels
=
4096
N_pixels
=
512
shape
=
[
N_pixels
]
*
dimensionality
signal_space
=
ift
.
RGSpace
(
shape
,
distances
=
L
/
N_pixels
)
if
use_nifty2go
:
harmonic_space
=
signal_space
.
get_default_codomain
()
else
:
harmonic_space
=
ift
.
FFTOperator
.
get_default_codomain
(
signal_space
)
harmonic_space
=
signal_space
.
get_default_codomain
()
fft
=
ift
.
FFTOperator
(
harmonic_space
,
target
=
signal_space
)
power_space
=
ift
.
PowerSpace
(
harmonic_space
)
...
...
@@ -56,12 +45,8 @@ if __name__ == "__main__":
power_spectrum
=
power_spectrum
)
np
.
random
.
seed
(
43
)
if
use_nifty2go
:
mock_power
=
ift
.
PS_field
(
power_space
,
power_spectrum
)
mock_harmonic
=
ift
.
power_synthesize
(
mock_power
,
real_signal
=
True
)
else
:
mock_power
=
ift
.
Field
(
power_space
,
val
=
power_spectrum
)
mock_harmonic
=
mock_power
.
power_synthesize
(
real_signal
=
True
)
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
)
...
...
@@ -69,19 +54,11 @@ if __name__ == "__main__":
R
=
ift
.
ResponseOperator
(
signal_space
,
sigma
=
(
response_sigma
,),
exposure
=
(
exposure
,))
data_domain
=
R
.
target
[
0
]
if
use_nifty2go
:
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
])
else
:
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
],
default_spaces
=
[
0
,
0
])
if
use_nifty2go
:
N
=
ift
.
DiagonalOperator
(
ift
.
Field
.
full
(
data_domain
,
mock_signal
.
var
()
/
signal_to_noise
).
weight
(
1
))
else
:
ndiag
=
ift
.
Field
(
data_domain
,
mock_signal
.
var
()
/
signal_to_noise
).
weight
(
1
)
N
=
ift
.
DiagonalOperator
(
data_domain
,
ndiag
)
R_harmonic
=
ift
.
ComposedOperator
([
fft
,
R
])
N
=
ift
.
DiagonalOperator
(
ift
.
Field
.
full
(
data_domain
,
mock_signal
.
var
()
/
signal_to_noise
).
weight
(
1
))
noise
=
ift
.
Field
.
from_random
(
domain
=
data_domain
,
random_type
=
'normal'
,
std
=
mock_signal
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
...
...
@@ -90,23 +67,12 @@ if __name__ == "__main__":
# Wiener filter
j
=
R_harmonic
.
adjoint_times
(
N
.
inverse_times
(
data
))
if
use_nifty2go
:
ctrl
=
ift
.
GradientNormController
(
verbose
=
True
,
iteration_limit
=
10
,
tol_abs_gradnorm
=
1e-4
/
nu
.
K
/
(
nu
.
m
**
(
0.5
*
dimensionality
)))
else
:
def
print_stats
(
a_energy
,
iteration
):
# returns current energy
x
=
a_energy
.
value
print
(
x
,
iteration
)
ctrl
=
ift
.
GradientNormController
(
callback
=
print_stats
,
iteration_limit
=
10
,
tol_abs_gradnorm
=
1e-4
/
nu
.
K
/
(
nu
.
m
**
(
0.5
*
dimensionality
)))
ctrl
=
ift
.
GradientNormController
(
verbose
=
True
,
tol_abs_gradnorm
=
1e-4
/
nu
.
K
/
(
nu
.
m
**
(
0.5
*
dimensionality
)))
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
inverter
=
ift
.
ConjugateGradient
(
controller
=
ctrl
)
if
use_nifty2go
:
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
)
wiener_curvature
=
ift
.
InversionEnabler
(
wiener_curvature
,
inverter
)
else
:
wiener_curvature
=
ift
.
library
.
WienerFilterCurvature
(
S
=
S
,
N
=
N
,
R
=
R_harmonic
,
inverter
=
inverter
)
wiener_curvature
=
ift
.
InversionEnabler
(
wiener_curvature
,
inverter
)
m
=
wiener_curvature
.
inverse_times
(
j
)
m_s
=
fft
(
m
)
...
...
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