Skip to content
GitLab
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
92422402
Commit
92422402
authored
Jul 12, 2017
by
Theo Steininger
Browse files
Merge branch 'working_on_demos' into 'master'
Working on demos See merge request !120
parents
c24c4f01
ac8f1fe2
Pipeline
#14786
passed with stages
in 12 minutes and 50 seconds
Changes
26
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
92422402
...
...
@@ -2,6 +2,11 @@
setup.cfg
.idea
.DS_Store
*.pyc
*.html
.document
.svn/
*.csv
# from https://github.com/github/gitignore/blob/master/Python.gitignore
...
...
@@ -95,4 +100,4 @@ ENV/
.spyderproject
# Rope project settings
.ropeproject
\ No newline at end of file
.ropeproject
demos/__init__.py
0 → 100644
View file @
92422402
demos/critical_filtering.py
0 → 100644
View file @
92422402
from
nifty
import
*
from
nifty.library.wiener_filter
import
WienerFilterEnergy
from
nifty.library.critical_filter
import
CriticalPowerEnergy
import
plotly.offline
as
pl
import
plotly.graph_objs
as
go
from
mpi4py
import
MPI
comm
=
MPI
.
COMM_WORLD
rank
=
comm
.
rank
np
.
random
.
seed
(
42
)
def
plot_parameters
(
m
,
t
,
p
,
p_d
):
x
=
log
(
t
.
domain
[
0
].
kindex
)
m
=
fft
.
adjoint_times
(
m
)
m
=
m
.
val
.
get_full_data
().
real
t
=
t
.
val
.
get_full_data
().
real
p
=
p
.
val
.
get_full_data
().
real
p_d
=
p_d
.
val
.
get_full_data
().
real
pl
.
plot
([
go
.
Heatmap
(
z
=
m
)],
filename
=
'map.html'
)
pl
.
plot
([
go
.
Scatter
(
x
=
x
,
y
=
t
),
go
.
Scatter
(
x
=
x
,
y
=
p
),
go
.
Scatter
(
x
=
x
,
y
=
p_d
)],
filename
=
"t.html"
)
class
AdjointFFTResponse
(
LinearOperator
):
def
__init__
(
self
,
FFT
,
R
,
default_spaces
=
None
):
super
(
AdjointFFTResponse
,
self
).
__init__
(
default_spaces
)
self
.
_domain
=
FFT
.
target
self
.
_target
=
R
.
target
self
.
R
=
R
self
.
FFT
=
FFT
def
_times
(
self
,
x
,
spaces
=
None
):
return
self
.
R
(
self
.
FFT
.
adjoint_times
(
x
))
def
_adjoint_times
(
self
,
x
,
spaces
=
None
):
return
self
.
FFT
(
self
.
R
.
adjoint_times
(
x
))
@
property
def
domain
(
self
):
return
self
.
_domain
@
property
def
target
(
self
):
return
self
.
_target
@
property
def
unitary
(
self
):
return
False
if
__name__
==
"__main__"
:
distribution_strategy
=
'not'
# Set up position space
s_space
=
RGSpace
([
128
,
128
])
# s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft
=
FFTOperator
(
s_space
)
h_space
=
fft
.
target
[
0
]
# Setting up power space
p_space
=
PowerSpace
(
h_space
,
logarithmic
=
True
,
distribution_strategy
=
distribution_strategy
)
# Choosing the prior correlation structure and defining correlation operator
p_spec
=
(
lambda
k
:
(.
5
/
(
k
+
1
)
**
3
))
S
=
create_power_operator
(
h_space
,
power_spectrum
=
p_spec
,
distribution_strategy
=
distribution_strategy
)
# Drawing a sample sh from the prior distribution in harmonic space
sp
=
Field
(
p_space
,
val
=
p_spec
,
distribution_strategy
=
distribution_strategy
)
sh
=
sp
.
power_synthesize
(
real_signal
=
True
)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument
=
DiagonalOperator
(
s_space
,
diagonal
=
1.
)
# Instrument._diagonal.val[200:400, 200:400] = 0
#Instrument._diagonal.val[64:512-64, 64:512-64] = 0
#Adding a harmonic transformation to the instrument
R
=
AdjointFFTResponse
(
fft
,
Instrument
)
noise
=
1.
N
=
DiagonalOperator
(
s_space
,
diagonal
=
noise
,
bare
=
True
)
n
=
Field
.
from_random
(
domain
=
s_space
,
random_type
=
'normal'
,
std
=
sqrt
(
noise
),
mean
=
0
)
# Creating the mock data
d
=
R
(
sh
)
+
n
# The information source
j
=
R
.
adjoint_times
(
N
.
inverse_times
(
d
))
realized_power
=
log
(
sh
.
power_analyze
(
logarithmic
=
p_space
.
config
[
"logarithmic"
],
nbin
=
p_space
.
config
[
"nbin"
]))
data_power
=
log
(
fft
(
d
).
power_analyze
(
logarithmic
=
p_space
.
config
[
"logarithmic"
],
nbin
=
p_space
.
config
[
"nbin"
]))
d_data
=
d
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
d_data
)],
filename
=
'data.html'
)
# minimization strategy
def
convergence_measure
(
a_energy
,
iteration
):
# returns current energy
x
=
a_energy
.
value
print
(
x
,
iteration
)
minimizer1
=
RelaxedNewton
(
convergence_tolerance
=
10e-2
,
convergence_level
=
2
,
iteration_limit
=
3
,
callback
=
convergence_measure
)
minimizer2
=
VL_BFGS
(
convergence_tolerance
=
0
,
iteration_limit
=
7
,
callback
=
convergence_measure
,
max_history_length
=
3
)
# Setting starting position
flat_power
=
Field
(
p_space
,
val
=
10e-8
)
m0
=
flat_power
.
power_synthesize
(
real_signal
=
True
)
t0
=
Field
(
p_space
,
val
=
log
(
1.
/
(
1
+
p_space
.
kindex
)
**
2
))
for
i
in
range
(
500
):
S0
=
create_power_operator
(
h_space
,
power_spectrum
=
exp
(
t0
),
distribution_strategy
=
distribution_strategy
)
# Initializing the nonlinear Wiener Filter energy
map_energy
=
WienerFilterEnergy
(
position
=
m0
,
d
=
d
,
R
=
R
,
N
=
N
,
S
=
S0
,
inverter
=
inverter
)
# Solving the Wiener Filter analytically
D0
=
map_energy
.
curvature
m0
=
D0
.
inverse_times
(
j
)
# Initializing the power energy with updated parameters
power_energy
=
CriticalPowerEnergy
(
position
=
t0
,
m
=
m0
,
D
=
D0
,
smoothness_prior
=
10.
,
samples
=
3
)
(
power_energy
,
convergence
)
=
minimizer1
(
power_energy
)
# Setting new power spectrum
t0
.
val
=
power_energy
.
position
.
val
.
real
# Plotting current estimate
print
i
if
i
%
50
==
0
:
plot_parameters
(
m0
,
t0
,
log
(
sp
),
data_power
)
demos/wiener_filter_advanced.py
0 → 100644
View file @
92422402
from
nifty
import
*
import
plotly.offline
as
pl
import
plotly.graph_objs
as
go
from
mpi4py
import
MPI
comm
=
MPI
.
COMM_WORLD
rank
=
comm
.
rank
np
.
random
.
seed
(
42
)
class
AdjointFFTResponse
(
LinearOperator
):
def
__init__
(
self
,
FFT
,
R
,
default_spaces
=
None
):
super
(
AdjointFFTResponse
,
self
).
__init__
(
default_spaces
)
self
.
_domain
=
FFT
.
target
self
.
_target
=
R
.
target
self
.
R
=
R
self
.
FFT
=
FFT
def
_times
(
self
,
x
,
spaces
=
None
):
return
self
.
R
(
self
.
FFT
.
adjoint_times
(
x
))
def
_adjoint_times
(
self
,
x
,
spaces
=
None
):
return
self
.
FFT
(
self
.
R
.
adjoint_times
(
x
))
@
property
def
domain
(
self
):
return
self
.
_domain
@
property
def
target
(
self
):
return
self
.
_target
@
property
def
unitary
(
self
):
return
False
if
__name__
==
"__main__"
:
distribution_strategy
=
'not'
# Set up position space
s_space
=
RGSpace
([
128
,
128
])
# s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft
=
FFTOperator
(
s_space
)
h_space
=
fft
.
target
[
0
]
# Setting up power space
p_space
=
PowerSpace
(
h_space
,
distribution_strategy
=
distribution_strategy
)
# Choosing the prior correlation structure and defining correlation operator
p_spec
=
(
lambda
k
:
(
42
/
(
k
+
1
)
**
3
))
S
=
create_power_operator
(
h_space
,
power_spectrum
=
p_spec
,
distribution_strategy
=
distribution_strategy
)
# Drawing a sample sh from the prior distribution in harmonic space
sp
=
Field
(
p_space
,
val
=
p_spec
,
distribution_strategy
=
distribution_strategy
)
sh
=
sp
.
power_synthesize
(
real_signal
=
True
)
ss
=
fft
.
adjoint_times
(
sh
)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05)
Instrument
=
DiagonalOperator
(
s_space
,
diagonal
=
1.
)
# Instrument._diagonal.val[200:400, 200:400] = 0
#Adding a harmonic transformation to the instrument
R
=
AdjointFFTResponse
(
fft
,
Instrument
)
signal_to_noise
=
1.
N
=
DiagonalOperator
(
s_space
,
diagonal
=
ss
.
var
()
/
signal_to_noise
,
bare
=
True
)
n
=
Field
.
from_random
(
domain
=
s_space
,
random_type
=
'normal'
,
std
=
ss
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
# Creating the mock data
d
=
R
(
sh
)
+
n
j
=
R
.
adjoint_times
(
N
.
inverse_times
(
d
))
# Choosing the minimization strategy
def
convergence_measure
(
energy
,
iteration
):
# returns current energy
x
=
energy
.
value
print
(
x
,
iteration
)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure)
minimizer
=
RelaxedNewton
(
convergence_tolerance
=
0
,
iteration_limit
=
1
,
callback
=
convergence_measure
)
#
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure,
# max_history_length=3)
#
inverter
=
ConjugateGradient
(
convergence_level
=
3
,
convergence_tolerance
=
10e-5
,
preconditioner
=
None
)
# Setting starting position
m0
=
Field
(
h_space
,
val
=
.
0
)
# Initializing the Wiener Filter energy
energy
=
WienerFilterEnergy
(
position
=
m0
,
d
=
d
,
R
=
R
,
N
=
N
,
S
=
S
,
inverter
=
inverter
)
D0
=
energy
.
curvature
# Solving the problem analytically
m0
=
D0
.
inverse_times
(
j
)
sample_variance
=
Field
(
sh
.
domain
,
val
=
0.
+
0j
)
sample_mean
=
Field
(
sh
.
domain
,
val
=
0.
+
0j
)
# sampling the uncertainty map
n_samples
=
1
for
i
in
range
(
n_samples
):
sample
=
sugar
.
generate_posterior_sample
(
m0
,
D0
)
sample_variance
+=
sample
**
2
sample_mean
+=
sample
variance
=
sample_variance
/
n_samples
-
(
sample_mean
/
n_samples
)
demos/wiener_filter.py
→
demos/wiener_filter
_easy
.py
View file @
92422402
...
...
@@ -19,7 +19,7 @@ if __name__ == "__main__":
correlation_length
=
0.1
#variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance
=
2.
#smoothing length
that
response (in same unit as L)
#smoothing length
of
response (in same unit as L)
response_sigma
=
0.1
#defining resolution (pixels per dimension)
...
...
@@ -62,15 +62,8 @@ if __name__ == "__main__":
d
=
R
(
ss
)
+
n
# Wiener filter
j
=
R
.
adjoint_times
(
N
.
inverse_times
(
d
))
D
=
PropagatorOperator
(
S
=
S
,
N
=
N
,
R
=
R
)
m
=
D
(
j
)
d_data
=
d
.
val
.
get_full_data
().
real
m_data
=
m
.
val
.
get_full_data
().
real
ss_data
=
ss
.
val
.
get_full_data
().
real
# if rank == 0:
# pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# pl.plot([go.Heatmap(z=m_data)], filename='map.html')
# pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')
demos/wiener_filter_hamiltonian.py
deleted
100644 → 0
View file @
c24c4f01
from
nifty
import
*
import
plotly.offline
as
pl
import
plotly.graph_objs
as
go
from
mpi4py
import
MPI
comm
=
MPI
.
COMM_WORLD
rank
=
comm
.
rank
np
.
random
.
seed
(
42
)
class
WienerFilterEnergy
(
Energy
):
def
__init__
(
self
,
position
,
D
,
j
):
# in principle not necessary, but useful in order to make the signature
# explicit
super
(
WienerFilterEnergy
,
self
).
__init__
(
position
)
self
.
D
=
D
self
.
j
=
j
def
at
(
self
,
position
):
return
self
.
__class__
(
position
,
D
=
self
.
D
,
j
=
self
.
j
)
@
property
def
value
(
self
):
D_inv_x
=
self
.
D_inverse_x
()
H
=
0.5
*
D_inv_x
.
vdot
(
self
.
position
)
-
self
.
j
.
vdot
(
self
.
position
)
return
H
.
real
@
property
def
gradient
(
self
):
D_inv_x
=
self
.
D_inverse_x
()
g
=
D_inv_x
-
self
.
j
return_g
=
g
.
copy_empty
(
dtype
=
np
.
float
)
return_g
.
val
=
g
.
val
.
real
return
return_g
@
property
def
curvature
(
self
):
class
Dummy
(
object
):
def
__init__
(
self
,
x
):
self
.
x
=
x
def
inverse_times
(
self
,
*
args
,
**
kwargs
):
return
self
.
x
.
times
(
*
args
,
**
kwargs
)
my_dummy
=
Dummy
(
self
.
D
)
return
my_dummy
@
memo
def
D_inverse_x
(
self
):
return
D
.
inverse_times
(
self
.
position
)
if
__name__
==
"__main__"
:
distribution_strategy
=
'not'
# Set up spaces and fft transformation
s_space
=
RGSpace
([
512
,
512
])
fft
=
FFTOperator
(
s_space
)
h_space
=
fft
.
target
[
0
]
p_space
=
PowerSpace
(
h_space
,
distribution_strategy
=
distribution_strategy
)
# create the field instances and power operator
pow_spec
=
(
lambda
k
:
(
42
/
(
k
+
1
)
**
3
))
S
=
create_power_operator
(
h_space
,
power_spectrum
=
pow_spec
,
distribution_strategy
=
distribution_strategy
)
sp
=
Field
(
p_space
,
val
=
lambda
z
:
pow_spec
(
z
)
**
(
1.
/
2
),
distribution_strategy
=
distribution_strategy
)
sh
=
sp
.
power_synthesize
(
real_signal
=
True
)
ss
=
fft
.
inverse_times
(
sh
)
# model the measurement process
R
=
SmoothingOperator
(
s_space
,
sigma
=
0.01
)
# R = DiagonalOperator(s_space, diagonal=1.)
# R._diagonal.val[200:400, 200:400] = 0
signal_to_noise
=
1
N
=
DiagonalOperator
(
s_space
,
diagonal
=
ss
.
var
()
/
signal_to_noise
,
bare
=
True
)
n
=
Field
.
from_random
(
domain
=
s_space
,
random_type
=
'normal'
,
std
=
ss
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
# create mock data
d
=
R
(
ss
)
+
n
# set up reconstruction objects
j
=
R
.
adjoint_times
(
N
.
inverse_times
(
d
))
D
=
PropagatorOperator
(
S
=
S
,
N
=
N
,
R
=
R
)
def
distance_measure
(
energy
,
iteration
):
x
=
energy
.
position
print
(
iteration
,
((
x
-
ss
).
norm
()
/
ss
.
norm
()).
real
)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=distance_measure)
minimizer
=
RelaxedNewton
(
convergence_tolerance
=
0
,
iteration_limit
=
2
,
callback
=
distance_measure
)
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# callback=distance_measure,
# max_history_length=3)
m0
=
Field
(
s_space
,
val
=
1.
)
energy
=
WienerFilterEnergy
(
position
=
m0
,
D
=
D
,
j
=
j
)
(
energy
,
convergence
)
=
minimizer
(
energy
)
m
=
energy
.
position
d_data
=
d
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
d_data
)],
filename
=
'data.html'
)
ss_data
=
ss
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
ss_data
)],
filename
=
'ss.html'
)
sh_data
=
sh
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
sh_data
)],
filename
=
'sh.html'
)
j_data
=
j
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
j_data
)],
filename
=
'j.html'
)
jabs_data
=
np
.
abs
(
j
.
val
.
get_full_data
())
jphase_data
=
np
.
angle
(
j
.
val
.
get_full_data
())
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
jabs_data
)],
filename
=
'j_abs.html'
)
pl
.
plot
([
go
.
Heatmap
(
z
=
jphase_data
)],
filename
=
'j_phase.html'
)
m_data
=
m
.
val
.
get_full_data
().
real
if
rank
==
0
:
pl
.
plot
([
go
.
Heatmap
(
z
=
m_data
)],
filename
=
'map.html'
)
# grad_data = grad.val.get_full_data().real
# if rank == 0:
# pl.plot([go.Heatmap(z=grad_data)], filename='grad.html')
demos/wiener_filter_harmonic.py
deleted
100644 → 0
View file @
c24c4f01
from
nifty
import
*
from
mpi4py
import
MPI
import
plotly.offline
as
py
import
plotly.graph_objs
as
go
comm
=
MPI
.
COMM_WORLD
rank
=
comm
.
rank
def
plot_maps
(
x
,
name
):
trace
=
[
None
]
*
len
(
x
)
keys
=
x
.
keys
()
field
=
x
[
keys
[
0
]]
domain
=
field
.
domain
[
0
]
shape
=
len
(
domain
.
shape
)
max_n
=
domain
.
shape
[
0
]
*
domain
.
distances
[
0
]
step
=
domain
.
distances
[
0
]
x_axis
=
np
.
arange
(
0
,
max_n
,
step
)
if
shape
==
1
:
for
ii
in
xrange
(
len
(
x
)):
trace
[
ii
]
=
go
.
Scatter
(
x
=
x_axis
,
y
=
x
[
keys
[
ii
]].
val
.
get_full_data
(),
name
=
keys
[
ii
])
fig
=
go
.
Figure
(
data
=
trace
)
py
.
plot
(
fig
,
filename
=
name
)
elif
shape
==
2
:
for
ii
in
xrange
(
len
(
x
)):
py
.
plot
([
go
.
Heatmap
(
z
=
x
[
keys
[
ii
]].
val
.
get_full_data
())],
filename
=
keys
[
ii
])
else
:
raise
TypeError
(
"Only 1D and 2D field plots are supported"
)
def
plot_power
(
x
,
name
):
layout
=
go
.
Layout
(
xaxis
=
dict
(
type
=
'log'
,
autorange
=
True
),
yaxis
=
dict
(
type
=
'log'
,
autorange
=
True
)
)
trace
=
[
None
]
*
len
(
x
)
keys
=
x
.
keys
()
field
=
x
[
keys
[
0
]]
domain
=
field
.
domain
[
0
]
x_axis
=
domain
.
kindex
for
ii
in
xrange
(
len
(
x
)):
trace
[
ii
]
=
go
.
Scatter
(
x
=
x_axis
,
y
=
x
[
keys
[
ii
]].
val
.
get_full_data
(),
name
=
keys
[
ii
])
fig
=
go
.
Figure
(
data
=
trace
,
layout
=
layout
)
py
.
plot
(
fig
,
filename
=
name
)
np
.
random
.
seed
(
42
)
if
__name__
==
"__main__"
:
distribution_strategy
=
'not'
# setting spaces
npix
=
np
.
array
([
500
])
# number of pixels
total_volume
=
1.
# total length
# setting signal parameters
lambda_s
=
.
05
# signal correlation length
sigma_s
=
10.
# signal variance
#setting response operator parameters
length_convolution
=
.
025
exposure
=
1.
# calculating parameters
k_0
=
4.
/
(
2
*
np
.
pi
*
lambda_s
)
a_s
=
sigma_s
**
2.
*
lambda_s
*
total_volume
# creation of spaces
# x1 = RGSpace([npix,npix], distances=total_volume / npix,
# zerocenter=False)
# k1 = RGRGTransformation.get_codomain(x1)
x1
=
HPSpace
(
64
)
k1
=
HPLMTransformation
.
get_codomain
(
x1
)
p1
=
PowerSpace
(
harmonic_partner
=
k1
,
logarithmic
=
False
)
# creating Power Operator with given spectrum
spec
=
(
lambda
k
:
a_s
/
(
1
+
(
k
/
k_0
)
**
2
)
**
2
)
p_field
=
Field
(
p1
,
val
=
spec
)
S_op
=
create_power_operator
(
k1
,
spec
)
# creating FFT-Operator and Response-Operator with Gaussian convolution
# adjust dtype_target probperly
Fft_op
=
FFTOperator
(
domain
=
x1
,
target
=
k1
,
domain_dtype
=
np
.
float64
,
target_dtype
=
np
.
float64
)
R_op
=
ResponseOperator
(
x1
,
sigma
=
[
length_convolution
],
exposure
=
[
exposure
])
# drawing a random field
sk
=
p_field
.
power_synthesize
(
real_power
=
True
,
mean
=
0.
)
s
=
Fft_op
.
adjoint_times
(
sk
)
signal_to_noise
=
1
N_op
=
DiagonalOperator
(
R_op
.
target
,
diagonal
=
s
.
var
()
/
signal_to_noise
,
bare
=
True
)
n
=
Field
.
from_random
(
domain
=
R_op
.
target
,
random_type
=
'normal'
,
std
=
s
.
std
()
/
np
.
sqrt
(
signal_to_noise
),
mean
=
0
)
d
=
R_op
(
s
)
+
n
# Wiener filter