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
609ba005
Commit
609ba005
authored
Apr 13, 2021
by
Martin Reinecke
Browse files
Merge branch 'nifty7_new_zm_handling' into 'NIFTy_7'
Revamp zero-mode handling See merge request
!607
parents
77d9d3e0
be037fd2
Pipeline
#100796
passed with stages
in 25 minutes and 56 seconds
Changes
8
Pipelines
22
Hide whitespace changes
Inline
Side-by-side
.gitlab-ci.yml
View file @
609ba005
...
...
@@ -110,6 +110,14 @@ run_getting_started_mf:
paths
:
-
'
*.png'
run_getting_density
:
stage
:
demo_runs
script
:
-
python3 demos/getting_started_density.py
artifacts
:
paths
:
-
'
*.png'
run_bernoulli
:
stage
:
demo_runs
script
:
...
...
ChangeLog.md
View file @
609ba005
Changes since NIFTy 6
=====================
New parametric amplitude model
------------------------------
The
`ift.CorrelatedFieldMaker`
now features two amplitude models. In addition
to the non-parametric one, one may choose to use a Matern kernel instead. The
method is aptly named
`add_fluctuations_matern`
. The major advantage of the
parametric model is its more intuitive scaling with the size of the position
space.
CorrelatedFieldMaker interface change
-------------------------------------
The interface of
`ift.CorrelatedFieldMaker.make`
and
`ift.CorrelatedFieldMaker.add_fluctuations`
changed; it now expects the mean
and the standard deviation of their various parameters not as separate
arguments but as a tuple.
The interface of
`ift.CorrelatedFieldMaker`
changed and instances of it may now
be instantiated directly without the previously required
`make`
method. Upon
initialization, no zero-mode must be specified as the normalization for the
different axes of the power respectively amplitude spectrum now only happens
once in the
`finalize`
method. There is now a new call named
`set_amplitude_total_offset`
to set the zero-mode. The method accepts either an
instance of
`ift.Operator`
or a tuple parameterizing a log-normal parameter.
Methods which require the zero-mode to be set raise a
`NotImplementedError`
if
invoked prior to having specified a zero-mode.
Furthermore, the interface of
`ift.CorrelatedFieldMaker.add_fluctuations`
changed; it now expects the mean and the standard deviation of their various
parameters not as separate arguments but as a tuple. The same applies to all
new and renamed methods of the
`CorrelatedFieldMaker`
class.
Furthermore, it is now possible to disable the asperity and the flexibility
together with the asperity in the correlated field model. Note that disabling
...
...
demos/getting_started_4_CorrelatedFields.ipynb
View file @
609ba005
...
...
@@ -87,7 +87,8 @@
"main_sample = None\n",
"def init_model(m_pars, fl_pars):\n",
" global main_sample\n",
" cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
" cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
" cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
" cf.add_fluctuations(**fl_pars)\n",
" field = cf.finalize(prior_info=0)\n",
" main_sample = ift.from_random(field.domain)\n",
...
...
@@ -95,7 +96,8 @@
" \n",
"# Helper: field and spectrum from parameter dictionaries + plotting\n",
"def eval_model(m_pars, fl_pars, title=None, samples=None):\n",
" cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
" cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
" cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
" cf.add_fluctuations(**fl_pars)\n",
" field = cf.finalize(prior_info=0)\n",
" spectrum = cf.amplitude\n",
...
...
@@ -473,7 +475,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.
8.3
"
"version": "3.
9.2
"
}
},
"nbformat": 4,
...
...
%% Cell type:markdown id: tags:
# Notebook showcasing the NIFTy 6 Correlated Field model
**Skip to `Parameter Showcases` for the meat/veggies ;)**
The field model roughly works like this:
`f = HT( A * zero_mode * xi ) + offset`
`A`
is a spectral power field which is constructed from power spectra that hold on subdomains of the target domain.
It is scaled by a zero mode operator and then pointwise multiplied by a gaussian excitation field, yielding
a representation of the field in harmonic space.
It is then transformed into the target real space and a offset added.
The power spectra
`A`
is constructed of are in turn constructed as the sum of a power law component
and an integrated Wiener process whose amplitude and roughness can be set.
## Setup code
%% Cell type:code id: tags:
```
python
import
nifty7
as
ift
import
matplotlib.pyplot
as
plt
import
numpy
as
np
ift
.
random
.
push_sseq_from_seed
(
43
)
n_pix
=
256
x_space
=
ift
.
RGSpace
(
n_pix
)
```
%% Cell type:code id: tags:
```
python
# Plotting routine
def
plot
(
fields
,
spectra
,
title
=
None
):
# Plotting preparation is normally handled by nifty7.Plot
# It is done manually here to be able to tweak details
# Fields are assumed to have identical domains
fig
=
plt
.
figure
(
tight_layout
=
True
,
figsize
=
(
12
,
3.5
))
if
title
is
not
None
:
fig
.
suptitle
(
title
,
fontsize
=
14
)
# Field
ax1
=
fig
.
add_subplot
(
1
,
2
,
1
)
ax1
.
axhline
(
y
=
0.
,
color
=
'k'
,
linestyle
=
'--'
,
alpha
=
0.25
)
for
field
in
fields
:
dom
=
field
.
domain
[
0
]
xcoord
=
np
.
arange
(
dom
.
shape
[
0
])
*
dom
.
distances
[
0
]
ax1
.
plot
(
xcoord
,
field
.
val
)
ax1
.
set_xlim
(
xcoord
[
0
],
xcoord
[
-
1
])
ax1
.
set_ylim
(
-
5.
,
5.
)
ax1
.
set_xlabel
(
'x'
)
ax1
.
set_ylabel
(
'f(x)'
)
ax1
.
set_title
(
'Field realizations'
)
# Spectrum
ax2
=
fig
.
add_subplot
(
1
,
2
,
2
)
for
spectrum
in
spectra
:
xcoord
=
spectrum
.
domain
[
0
].
k_lengths
ycoord
=
spectrum
.
val_rw
()
ycoord
[
0
]
=
ycoord
[
1
]
ax2
.
plot
(
xcoord
,
ycoord
)
ax2
.
set_ylim
(
1e-6
,
10.
)
ax2
.
set_xscale
(
'log'
)
ax2
.
set_yscale
(
'log'
)
ax2
.
set_xlabel
(
'k'
)
ax2
.
set_ylabel
(
'p(k)'
)
ax2
.
set_title
(
'Power Spectrum'
)
fig
.
align_labels
()
plt
.
show
()
# Helper: draw main sample
main_sample
=
None
def
init_model
(
m_pars
,
fl_pars
):
global
main_sample
cf
=
ift
.
CorrelatedFieldMaker
.
make
(
**
m_pars
)
cf
=
ift
.
CorrelatedFieldMaker
(
m_pars
[
"prefix"
])
cf
.
set_amplitude_total_offset
(
m_pars
[
"offset_mean"
],
m_pars
[
"offset_std"
])
cf
.
add_fluctuations
(
**
fl_pars
)
field
=
cf
.
finalize
(
prior_info
=
0
)
main_sample
=
ift
.
from_random
(
field
.
domain
)
print
(
"model domain keys:"
,
field
.
domain
.
keys
())
# Helper: field and spectrum from parameter dictionaries + plotting
def
eval_model
(
m_pars
,
fl_pars
,
title
=
None
,
samples
=
None
):
cf
=
ift
.
CorrelatedFieldMaker
.
make
(
**
m_pars
)
cf
=
ift
.
CorrelatedFieldMaker
(
m_pars
[
"prefix"
])
cf
.
set_amplitude_total_offset
(
m_pars
[
"offset_mean"
],
m_pars
[
"offset_std"
])
cf
.
add_fluctuations
(
**
fl_pars
)
field
=
cf
.
finalize
(
prior_info
=
0
)
spectrum
=
cf
.
amplitude
if
samples
is
None
:
samples
=
[
main_sample
]
field_realizations
=
[
field
(
s
)
for
s
in
samples
]
spectrum_realizations
=
[
spectrum
.
force
(
s
)
for
s
in
samples
]
plot
(
field_realizations
,
spectrum_realizations
,
title
)
def
gen_samples
(
key_to_vary
):
if
key_to_vary
is
None
:
return
[
main_sample
]
dct
=
main_sample
.
to_dict
()
subdom_to_vary
=
dct
.
pop
(
key_to_vary
).
domain
samples
=
[]
for
i
in
range
(
8
):
d
=
dct
.
copy
()
d
[
key_to_vary
]
=
ift
.
from_random
(
subdom_to_vary
)
samples
.
append
(
ift
.
MultiField
.
from_dict
(
d
))
return
samples
def
vary_parameter
(
parameter_key
,
values
,
samples_vary_in
=
None
):
s
=
gen_samples
(
samples_vary_in
)
for
v
in
values
:
if
parameter_key
in
cf_make_pars
.
keys
():
m_pars
=
{
**
cf_make_pars
,
parameter_key
:
v
}
eval_model
(
m_pars
,
cf_x_fluct_pars
,
f
"
{
parameter_key
}
=
{
v
}
"
,
s
)
else
:
fl_pars
=
{
**
cf_x_fluct_pars
,
parameter_key
:
v
}
eval_model
(
cf_make_pars
,
fl_pars
,
f
"
{
parameter_key
}
=
{
v
}
"
,
s
)
```
%% Cell type:markdown id: tags:
## Before the Action: The Moment-Matched Log-Normal Distribution
Many properties of the correlated field are modelled as being lognormally distributed.
The distribution models are parametrized via their means and standard-deviations (first and second position in tuple).
To get a feeling of how the ratio of the
`mean`
and
`stddev`
parameters influences the distribution shape,
here are a few example histograms: (observe the x-axis!)
%% Cell type:code id: tags:
```
python
fig
=
plt
.
figure
(
figsize
=
(
13
,
3.5
))
mean
=
1.0
sigmas
=
[
1.0
,
0.5
,
0.1
]
for
i
in
range
(
3
):
op
=
ift
.
LognormalTransform
(
mean
=
mean
,
sigma
=
sigmas
[
i
],
key
=
'foo'
,
N_copies
=
0
)
op_samples
=
np
.
array
(
[
op
(
s
).
val
for
s
in
[
ift
.
from_random
(
op
.
domain
)
for
i
in
range
(
10000
)]])
ax
=
fig
.
add_subplot
(
1
,
3
,
i
+
1
)
ax
.
hist
(
op_samples
,
bins
=
50
)
ax
.
set_title
(
f
"mean =
{
mean
}
, sigma =
{
sigmas
[
i
]
}
"
)
ax
.
set_xlabel
(
'x'
)
del
op_samples
plt
.
show
()
```
%% Cell type:markdown id: tags:
## The Neutral Field
To demonstrate the effect of all parameters, first a 'neutral' set of parameters
is defined which then are varied one by one, showing the effect of the variation
on the generated field realizations and the underlying power spectrum from which
they were drawn.
As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen.
%% Cell type:code id: tags:
```
python
# Neutral model parameters yielding a quasi-constant field
cf_make_pars
=
{
'offset_mean'
:
0.
,
'offset_std'
:
(
1e-3
,
1e-16
),
'prefix'
:
''
}
cf_x_fluct_pars
=
{
'target_subdomain'
:
x_space
,
'fluctuations'
:
(
1e-3
,
1e-16
),
'flexibility'
:
(
1e-3
,
1e-16
),
'asperity'
:
(
1e-3
,
1e-16
),
'loglogavgslope'
:
(
0.
,
1e-16
)
}
init_model
(
cf_make_pars
,
cf_x_fluct_pars
)
```
%% Cell type:code id: tags:
```
python
# Show neutral field
eval_model
(
cf_make_pars
,
cf_x_fluct_pars
,
"Neutral Field"
)
```
%% Cell type:markdown id: tags:
# Parameter Showcases
## The `fluctuations` parameters of `add_fluctuations()`
determine the
**amplitude of variations along the field dimension**
for which
`add_fluctuations`
is called.
`fluctuations[0]`
set the _average_ amplitude of the fields fluctuations along the given dimension,
\
`fluctuations[1]`
sets the width and shape of the amplitude distribution.
The amplitude is modelled as being log-normally distributed,
see
`The Moment-Matched Log-Normal Distribution`
above for details.
#### `fluctuations` mean:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'fluctuations'
,
[(
0.05
,
1e-16
),
(
0.5
,
1e-16
),
(
2.
,
1e-16
)],
samples_vary_in
=
'xi'
)
```
%% Cell type:markdown id: tags:
#### `fluctuations` std:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'fluctuations'
,
[(
1.
,
0.01
),
(
1.
,
0.1
),
(
1.
,
1.
)],
samples_vary_in
=
'fluctuations'
)
cf_x_fluct_pars
[
'fluctuations'
]
=
(
1.
,
1e-16
)
```
%% Cell type:markdown id: tags:
## The `loglogavgslope` parameters of `add_fluctuations()`
determine
**the slope of the loglog-linear (power law) component of the power spectrum**
.
The slope is modelled to be normally distributed.
#### `loglogavgslope` mean:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'loglogavgslope'
,
[(
-
6.
,
1e-16
),
(
-
2.
,
1e-16
),
(
2.
,
1e-16
)],
samples_vary_in
=
'xi'
)
```
%% Cell type:markdown id: tags:
#### `loglogavgslope` std:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'loglogavgslope'
,
[(
-
2.
,
0.02
),
(
-
2.
,
0.2
),
(
-
2.
,
2.0
)],
samples_vary_in
=
'loglogavgslope'
)
cf_x_fluct_pars
[
'loglogavgslope'
]
=
(
-
2.
,
1e-16
)
```
%% Cell type:markdown id: tags:
## The `flexibility` parameters of `add_fluctuations()`
determine
**the amplitude of the integrated Wiener process component of the power spectrum**
(how strong the power spectrum varies besides the power-law).
`flexibility[0]`
sets the _average_ amplitude of the i.g.p. component,
\
`flexibility[1]`
sets how much the amplitude can vary.
\
These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior.
#### `flexibility` mean:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'flexibility'
,
[(
0.4
,
1e-16
),
(
4.0
,
1e-16
),
(
12.0
,
1e-16
)],
samples_vary_in
=
'spectrum'
)
```
%% Cell type:markdown id: tags:
#### `flexibility` std:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'flexibility'
,
[(
4.
,
0.02
),
(
4.
,
0.2
),
(
4.
,
2.
)],
samples_vary_in
=
'flexibility'
)
cf_x_fluct_pars
[
'flexibility'
]
=
(
4.
,
1e-16
)
```
%% Cell type:markdown id: tags:
## The `asperity` parameters of `add_fluctuations()`
`asperity`
determines
**how rough the integrated Wiener process component of the power spectrum is**
.
`asperity[0]`
sets the average roughness,
`asperity[1]`
sets how much the roughness can vary.
\
These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior.
#### `asperity` mean:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'asperity'
,
[(
0.001
,
1e-16
),
(
1.0
,
1e-16
),
(
5.
,
1e-16
)],
samples_vary_in
=
'spectrum'
)
```
%% Cell type:markdown id: tags:
#### `asperity` std:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'asperity'
,
[(
1.
,
0.01
),
(
1.
,
0.1
),
(
1.
,
1.
)],
samples_vary_in
=
'asperity'
)
cf_x_fluct_pars
[
'asperity'
]
=
(
1.
,
1e-16
)
```
%% Cell type:markdown id: tags:
## The `offset_mean` parameter of `CorrelatedFieldMaker.make()`
The
`offset_mean`
parameter defines a global additive offset on the field realizations.
If the field is used for a lognormal model
`f = field.exp()`
, this acts as a global signal magnitude offset.
%% Cell type:code id: tags:
```
python
# Reset model to neutral
cf_x_fluct_pars
[
'fluctuations'
]
=
(
1e-3
,
1e-16
)
cf_x_fluct_pars
[
'flexibility'
]
=
(
1e-3
,
1e-16
)
cf_x_fluct_pars
[
'asperity'
]
=
(
1e-3
,
1e-16
)
cf_x_fluct_pars
[
'loglogavgslope'
]
=
(
1e-3
,
1e-16
)
```
%% Cell type:code id: tags:
```
python
vary_parameter
(
'offset_mean'
,
[
3.
,
0.
,
-
2.
])
```
%% Cell type:markdown id: tags:
## The `offset_std` parameters of `CorrelatedFieldMaker.make()`
Variation of the global offset of the field are modelled as being log-normally distributed.
See
`The Moment-Matched Log-Normal Distribution`
above for details.
The
`offset_std[0]`
parameter sets how much NIFTy will vary the offset
*on average*
.
\
The
`offset_std[1]`
parameters defines the with and shape of the offset variation distribution.
#### `offset_std` mean:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'offset_std'
,
[(
1e-16
,
1e-16
),
(
0.5
,
1e-16
),
(
2.
,
1e-16
)],
samples_vary_in
=
'xi'
)
```
%% Cell type:markdown id: tags:
#### `offset_std` std:
%% Cell type:code id: tags:
```
python
vary_parameter
(
'offset_std'
,
[(
1.
,
0.01
),
(
1.
,
0.1
),
(
1.
,
1.
)],
samples_vary_in
=
'zeromode'
)
```
...
...
demos/getting_started_5_mf.py
View file @
609ba005
...
...
@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-202
0
Max-Planck-Society
# Copyright(C) 2013-202
1
Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...
...
@@ -73,14 +73,17 @@ def main():
sp2
=
ift
.
RGSpace
(
npix2
)
# Set up signal model
cfmaker
=
ift
.
CorrelatedFieldMaker
.
make
(
0.
,
(
1e-2
,
1e-6
),
''
)
cfmaker
.
add_fluctuations
(
sp1
,
(
0.1
,
1e-2
),
(
2
,
.
2
),
(.
01
,
.
5
),
(
-
4
,
2.
),
'amp1'
)
cfmaker
.
add_fluctuations
(
sp2
,
(
0.1
,
1e-2
),
(
2
,
.
2
),
(.
01
,
.
5
),
(
-
3
,
1
),
'amp2'
)
cfmaker
=
ift
.
CorrelatedFieldMaker
(
''
)
cfmaker
.
add_fluctuations
(
sp1
,
(
0.1
,
1e-2
),
(
2
,
.
2
),
(.
01
,
.
5
),
(
-
4
,
2.
),
'amp1'
)
cfmaker
.
add_fluctuations
(
sp2
,
(
0.1
,
1e-2
),
(
2
,
.
2
),
(.
01
,
.
5
),
(
-
3
,
1
),
'amp2'
)
cfmaker
.
set_amplitude_total_offset
(
0.
,
(
1e-2
,
1e-6
))
correlated_field
=
cfmaker
.
finalize
()
pspec1
=
cfmaker
.
normalized_amplitudes
[
0
]
**
2
pspec2
=
cfmaker
.
normalized_amplitudes
[
1
]
**
2
normalized_amp
=
cfmaker
.
get_normalized_amplitudes
()
pspec1
=
normalized_amp
[
0
]
**
2
pspec2
=
normalized_amp
[
1
]
**
2
DC
=
SingleDomain
(
correlated_field
.
target
,
position_space
)
# Apply a nonlinearity
...
...
demos/getting_started_density.py
0 → 100644
View file @
609ba005
#!/usr/bin/env python3
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2021 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
############################################################
# Density estimation
#
# Compute a density estimate for a log-normal process measured by a
# Poissonian likelihood.
#
# Demo takes a while to compute
#############################################################
import
numpy
as
np
import
nifty7
as
ift
def
density_estimator
(
domain
,
pad
=
1.0
,
cf_fluctuations
=
None
,
cf_azm_uniform
=
None
):
cf_azm_uniform_sane_default
=
(
1e-4
,
1.0
)
cf_fluctuations_sane_default
=
{
"scale"
:
(
0.5
,
0.3
),
"cutoff"
:
(
4.0
,
3.0
),
"loglogslope"
:
(
-
6.0
,
3.0
)
}
domain
=
ift
.
DomainTuple
.
make
(
domain
)
dom_scaling
=
1.
+
np
.
broadcast_to
(
pad
,
(
len
(
domain
.
axes
),
))
if
cf_fluctuations
is
None
:
cf_fluctuations
=
cf_fluctuations_sane_default
if
cf_azm_uniform
is
None
:
cf_azm_uni
=
cf_azm_uniform_sane_default
domain_padded
=
[]
for
d_scl
,
d
in
zip
(
dom_scaling
,
domain
):
if
not
isinstance
(
d
,
ift
.
RGSpace
)
or
d
.
harmonic
:
te
=
[
f
"unexpected domain encountered in `domain`:
{
domain
}
"
]
te
+=
"expected a non-harmonic `ift.RGSpace`"
raise
TypeError
(
"
\n
"
.
join
(
te
))
shape_padded
=
tuple
((
d_scl
*
np
.
array
(
d
.
shape
)).
astype
(
int
))
domain_padded
.
append
(
ift
.
RGSpace
(
shape_padded
,
distances
=
d
.
distances
))
domain_padded
=
ift
.
DomainTuple
.
make
(
domain_padded
)
# Set up the signal model
azm_offset_mean
=
0.0
# The zero-mode should be inferred only from the data
cfmaker
=
ift
.
CorrelatedFieldMaker
(
""
)
for
i
,
d
in
enumerate
(
domain_padded
):
if
isinstance
(
cf_fluctuations
,
(
list
,
tuple
)):
cf_fl
=
cf_fluctuations
[
i
]
else
:
cf_fl
=
cf_fluctuations
cfmaker
.
add_fluctuations_matern
(
d
,
**
cf_fl
,
prefix
=
f
"ax
{
i
}
"
)
scalar_domain
=
ift
.
DomainTuple
.
scalar_domain
()
uniform
=
ift
.
UniformOperator
(
scalar_domain
,
*
cf_azm_uni
)
azm
=
uniform
.
ducktape
(
"zeromode"
)
cfmaker
.
set_amplitude_total_offset
(
azm_offset_mean
,
azm
)
correlated_field
=
cfmaker
.
finalize
(
0
).
clip
(
-
10.
,
10.
)
normalized_amplitudes
=
cfmaker
.
get_normalized_amplitudes
()
domain_shape
=
tuple
(
d
.
shape
for
d
in
domain
)
slc
=
ift
.
SliceOperator
(
correlated_field
.
target
,
domain_shape
)
signal
=
ift
.
exp
(
slc
@
correlated_field
)
model_operators
=
{
"correlated_field"
:
correlated_field
,
"select_subset"
:
slc
,
"amplitude_total_offset"
:
azm
,
"normalized_amplitudes"
:
normalized_amplitudes
}
return
signal
,
model_operators
if
__name__
==
"__main__"
:
filename
=
"getting_started_density_{}.png"
ift
.
random
.
push_sseq_from_seed
(
42
)
# Set up signal domain
npix1
=
128
npix2
=
128
sp1
=
ift
.
RGSpace
(
npix1
)
sp2
=
ift
.
RGSpace
(
npix2
)
position_space
=
ift
.
DomainTuple
.
make
((
sp1
,
sp2
))
signal
,
ops
=
density_estimator
(
position_space
)
correlated_field
=
ops
[
"correlated_field"
]
data_space
=
signal
.
target
# Generate mock signal and data
rng
=
ift
.
random
.
current_rng
()
mock_position
=
ift
.
from_random
(
signal
.
domain
,
"normal"
)
data
=
ift
.
Field
.
from_raw
(
data_space
,
rng
.
poisson
(
signal
(
mock_position
).
val
))
# Rejoin domains for plotting
plotting_domain
=
ift
.
DomainTuple
.
make
(
ift
.
RGSpace
((
npix1
,
npix2
)))
plotting_domain_expanded
=
ift
.
DomainTuple
.
make
(
ift
.
RGSpace
((
2
*
npix1
,
2
*
npix2
)))
plot
=
ift
.
Plot
()
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain_expanded
,
ift
.
exp
(
correlated_field
(
mock_position
)).
val
),
title
=
"Pre-Slicing Truth"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
signal
(
mock_position
).
val
),
title
=
"Ground Truth"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
data
.
val
),
title
=
"Data"
)
plot
.
output
(
ny
=
1
,
nx
=
3
,
xsize
=
10
,
ysize
=
3
,
name
=
filename
.
format
(
"setup"
))
print
(
"Setup saved as"
,
filename
.
format
(
"setup"
))
# Minimization parameters
ic_sampling
=
ift
.
AbsDeltaEnergyController
(
name
=
"Sampling"
,
deltaE
=
0.01
,
iteration_limit
=
100
)
ic_newton
=
ift
.
AbsDeltaEnergyController
(
name
=
"Newton"
,
deltaE
=
0.01
,
iteration_limit
=
35
)
ic_sampling
.
enable_logging
()
ic_newton
.
enable_logging
()
minimizer
=
ift
.
NewtonCG
(
ic_newton
,
enable_logging
=
True
)
# Number of samples used to estimate the KL
n_samples
=
5
# Set up likelihood and information Hamiltonian
likelihood
=
ift
.
PoissonianEnergy
(
data
)
@
signal
ham
=
ift
.
StandardHamiltonian
(
likelihood
,
ic_sampling
)
# Start minimization
initial_mean
=
ift
.
MultiField
.
full
(
ham
.
domain
,
0.
)
mean
=
initial_mean
for
i
in
range
(
5
):
# Draw new samples and minimize KL
kl
=
ift
.
MetricGaussianKL
.
make
(
mean
,
ham
,
n_samples
,
True
)
kl
,
convergence
=
minimizer
(
kl
)
mean
=
kl
.
position
# Plot current reconstruction
plot
=
ift
.
Plot
()
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain_expanded
,
ift
.
exp
(
correlated_field
(
mock_position
)).
val
),
title
=
"Ground truth"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
signal
(
mock_position
).
val
),
title
=
"Ground truth"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
signal
(
kl
.
position
).
val
),
title
=
"Reconstruction"
,
)
plot
.
add
(
(
ic_newton
.
history
,
ic_sampling
.
history
,
minimizer
.
inversion_history
),
label
=
[
"kl"
,
"Sampling"
,
"Newton inversion"
],
title
=
"Cumulative energies"
,
s
=
[
None
,
None
,
1
],
alpha
=
[
None
,
0.2
,
None
],
)
plot
.
output
(
nx
=
3
,
ny
=
2
,
ysize
=
10
,
xsize
=
15
,
name
=
filename
.
format
(
f
"loop_
{
i
:
02
d
}
"
)
)
# Done, draw posterior samples
sc
=
ift
.
StatCalculator
()
sc_unsliced
=
ift
.
StatCalculator
()
for
sample
in
kl
.
samples
:
sc
.
add
(
signal
(
sample
+
kl
.
position
))
sc_unsliced
.
add
(
ift
.
exp
(
correlated_field
(
sample
+
kl
.
position
)))
# Plotting
plot
=
ift
.
Plot
()
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
sc
.
mean
.
val
),
title
=
"Posterior Mean"
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain
,
ift
.
sqrt
(
sc
.
var
).
val
),
title
=
"Posterior Standard Deviation"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain_expanded
,
sc_unsliced
.
mean
.
val
),
title
=
"Posterior Unsliced Mean"
,
)
plot
.
add
(
ift
.
Field
.
from_raw
(
plotting_domain_expanded
,
ift
.
sqrt
(
sc_unsliced
.
var
).
val
),
title
=
"Posterior Unsliced Standard Deviation"
,
)
plot
.
output
(
xsize
=
15
,
ysize
=
15
,
name
=
filename
.
format
(
"results"
))
print
(
"Saved results as"
,
filename
.
format
(
"results"
))
src/library/correlated_fields.py
View file @
609ba005
...
...
@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-202
0
Max-Planck-Society
# Copyright(C) 2013-202
1
Max-Planck-Society
# Authors: Philipp Frank, Philipp Arras, Philipp Haim;
# Matern Kernel by Matteo Guardiani, Jakob Roth and
# Gordian Edenhofer
...
...
@@ -37,6 +37,7 @@ from ..operators.distributors import PowerDistributor
from
..operators.endomorphic_operator
import
EndomorphicOperator
from
..operators.harmonic_operators
import
HarmonicTransformOperator
from
..operators.linear_operator
import
LinearOperator
from
..operators.mask_operator
import
MaskOperator
from
..operators.normal_operators
import
LognormalTransform
,
NormalTransform
from
..operators.operator
import
Operator
from
..operators.simple_linear_operators
import
VdotOperator
,
ducktape
...
...
@@ -188,12 +189,13 @@ class _Normalization(Operator):
def
apply
(
self
,
x
):
self
.
_check_input
(
x
)
spec
=
x
.
ptw
(
"exp"
)
# NOTE, this "normalizes" also the zero-mode which is supposed to be
# left untouched by this operator. This is not a problem since this
# zero mode is multiplied with zero later on in the model and modelled
# differently.
# NOTE, see the note in the doc-string on why this is not a proper
# normalization!
# NOTE, this "normalizes" also the zero-mode which is supposed to be
# left untouched by this operator. Since the multiplicity of the
# zero-mode is set to 0, the norm does not contain traces of it.
# However, it wrongly sets the zeroth entry of the result. Luckily,
# in subsequent calls, the zeroth entry is not used in the CF model.
return
(
self
.
_specsum
(
spec
).
reciprocal
()
*
spec
).
sqrt
()
...
...
@@ -227,7 +229,7 @@ class _Distributor(LinearOperator):
class
_AmplitudeMatern
(
Operator
):
def
__init__
(
self
,
pow_spc
,
scale
,
cutoff
,
loglogslope
,
azm
,
totvol
):
def
__init__
(
self
,
pow_spc
,
scale
,
cutoff
,
loglogslope
,
totvol
):
expander
=
ContractionOperator
(
pow_spc
,
spaces
=
None
).
adjoint
k_squared
=
makeField
(
pow_spc
,
pow_spc
.
k_lengths
**
2
)
...
...
@@ -253,7 +255,7 @@ class _AmplitudeMatern(Operator):
vol1
[
1
:]
=
totvol
**
0.5
vol0
=
Adder
(
makeField
(
pow_spc
,
vol0
))
vol1
=
DiagonalOperator
(
makeField
(
pow_spc
,
vol1
))
op
=
vol0
@
(
vol1
@
(
expander
@
azm
.
ptw
(
"reciprocal"
))
*
op
)
op
=
vol0
@
vol1
@
op
# std = sqrt of integral of power spectrum
self
.
_fluc
=
op
.
power
(
2
).
integrate
().
sqrt
()
...
...
@@ -271,7 +273,7 @@ class _AmplitudeMatern(Operator):
class
_Amplitude
(
Operator
):
def
__init__
(
self
,
target
,
fluctuations
,
flexibility
,
asperity
,
loglogavgslope
,
azm
,
totvol
,
key
,
dofdex
):
loglogavgslope
,
totvol
,
key
,
dofdex
):
"""
fluctuations > 0
flexibility > 0 or None
...
...
@@ -294,7 +296,6 @@ class _Amplitude(Operator):
N_copies
=
0
space
=
0
distributed_tgt
=
target
=
makeDomain
(
target
)
azm_expander
=
ContractionOperator
(
distributed_tgt
,
spaces
=
space
).
adjoint
assert
isinstance
(
target
[
space
],
PowerSpace
)
twolog
=
_TwoLogIntegrations
(
target
,
space
)
...
...
@@ -359,11 +360,11 @@ class _Amplitude(Operator):
if
N_copies
>
0
:
op
=
Distributor
@
op
sig_fluc
=
Distributor
@
sig_fluc
op
=
Adder
(
Distributor
(
vol0
))
@
(
sig_fluc
*
(
azm_expander
@
azm
.
ptw
(
"reciprocal"
))
*
op
)
op
=
Adder
(
Distributor
(
vol0
))
@
(
sig_fluc
*
op
)
self
.
_fluc
=
(
_Distributor
(
dofdex
,
fluctuations
.
target
,
distributed_tgt
[
0
])
@
fluctuations
)
else
:
op
=
Adder
(
vol0
)
@
(
sig_fluc
*
(
azm_expander
@
azm
.
ptw
(
"reciprocal"
))
*
op
)
op
=
Adder
(
vol0
)
@
(
sig_fluc
*
op
)
self
.
_fluc
=
fluctuations
self
.
apply
=
op
.
apply
...
...
@@ -384,11 +385,10 @@ class CorrelatedFieldMaker:
The correlated field models are parametrized by creating
square roots of power spectrum operators ("amplitudes") via calls to
:func:`add_fluctuations` that act on the targeted field subdomains.
During creation of the :class:`CorrelatedFieldMaker` via
:func:`make`, a global offset from zero of the field model
can be defined and an operator applying fluctuations
around this offset is parametrized.
:func:`add_fluctuations*` that act on the targeted field subdomains.
During creation of the :class:`CorrelatedFieldMaker`, a global
offset from zero of the field model can be defined and an operator
applying fluctuations around this offset is parametrized.
The resulting correlated field model operator has a
:class:`~nifty7.multi_domain.MultiDomain` as its domain and
...
...
@@ -403,39 +403,21 @@ class CorrelatedFieldMaker:
:func:`finalize`, which returns the configured operator.
An operator representing an array of correlated field models
can be constructed by setting the `total_N` parameter of
:func:`make`. It will have an
:class:`~nifty7.domains.unstructured_domain.UnstructuredDomain`
can be constructed by setting the `total_N` parameter of. It will
have an :class:`~nifty7.domains.unstructured_domain.UnstructuredDomain`
of shape `(total_N,)` prepended to its target domain and represent
`total_N` correlated fields simulataneously.
The degree of information sharing between the correlated field
models can be configured via the `dofdex` parameters
of :func:`make` and :func:`add_fluctuations`.
See the methods :func:`make`, :func:`add_fluctuations`
and :func:`finalize` for further usage information."""
def
__init__
(
self
,
offset_mean
,
offset_fluctuations_op
,
prefix
,
total_N
):
if
not
isinstance
(
offset_fluctuations_op
,
Operator
):
raise
TypeError
(
"offset_fluctuations_op needs to be an operator"
)
self
.
_a
=
[]
self
.
_target_subdomains
=
[]
models can be configured via the `dofdex` parameter of
:func:`add_fluctuations`.
self
.
_offset_mean
=
offset_mean
self
.
_azm
=
offset_fluctuations_op
self
.
_prefix
=
prefix
self
.
_total_N
=
total_N
@
staticmethod
def
make
(
offset_mean
,
offset_std
,
prefix
,
total_N
=
0
,
dofdex
=
None
):
"""Returns a CorrelatedFieldMaker object.
See the methods :func:`add_fluctuations*` and :func:`finalize` for
further usage information."""
def
__init__
(
self
,
prefix
,
total_N
=
0
):
"""Instantiate a CorrelatedFieldMaker object.
Parameters
----------
offset_mean : float
Mean offset from zero of the correlated field to be made.
offset_std : tuple of float
Mean standard deviation and standard deviation of the standard
deviation of the offset. No, this is not a word duplication.
prefix : string
Prefix to the names of the domains of the cf operator to be made.
This determines the names of the operator domain.
...
...
@@ -444,26 +426,14 @@ class CorrelatedFieldMaker:
If not 0, the first entry of the operators target will be an
:class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
with length `total_N`.
dofdex : np.array of integers, optional
An integer array specifying the zero mode models used if
total_N > 1. It needs to have length of total_N. If total_N=3 and
dofdex=[0,0,1], that means that two models for the zero mode are
instantiated; the first one is used for the first and second
field model and the second is used for the third field model.
*If not specified*, use the same zero mode model for all
constructed field models.
"""
if
dofdex
is
None
:
dofdex
=
np
.
full
(
total_N
,
0
)
elif
len
(
dofdex
)
!=
total_N
:
raise
ValueError
(
"length of dofdex needs to match total_N"
)
N
=
max
(
dofdex
)
+
1
if
total_N
>
0
else
0
if
len
(
offset_std
)
!=
2
:
raise
TypeError
zm
=
LognormalTransform
(
*
offset_std
,
prefix
+
'zeromode'
,
N
)
if
total_N
>
0
:
zm
=
_Distributor
(
dofdex
,
zm
.
target
,
UnstructuredDomain
(
total_N
))
@
zm
return
CorrelatedFieldMaker
(
offset_mean
,
zm
,
prefix
,
total_N
)
self
.
_azm
=
None
self
.
_offset_mean
=
None
self
.
_a
=
[]
self
.
_target_subdomains
=
[]
self
.
_prefix
=
prefix
self
.
_total_N
=
total_N
def
add_fluctuations
(
self
,
target_subdomain
,
...
...
@@ -573,7 +543,7 @@ class CorrelatedFieldMaker:
avgsl
=
NormalTransform
(
*
loglogavgslope
,
pre
+
'loglogavgslope'
,
N
)