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
04c52785
Commit
04c52785
authored
Jun 19, 2020
by
Philipp Arras
Browse files
Merge branch 'nifty627' into 'NIFTy_7'
Nifty627 See merge request
!543
parents
de268998
0ba420fb
Pipeline
#76997
passed with stages
in 12 minutes and 33 seconds
Changes
8
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
ChangeLog.md
View file @
04c52785
...
...
@@ -76,6 +76,19 @@ print(met)
print(met.draw_sample())
```
New approach for sampling complex numbers
=========================================
When calling draw_sample_with_dtype with a complex dtype,
the variance is now used for the imaginary part and real part separately.
This is done in order to be consistent with the Hamiltonian.
Note that by this,
```
np.std(ift.from_random(domain, 'normal', dtype=np.complex128).val)
```
`
does not give 1, but sqrt(2) as a result.
MPI parallelisation over samples in MetricGaussianKL
----------------------------------------------------
...
...
@@ -91,6 +104,7 @@ the generation of reproducible random numbers in the presence of MPI parallelism
and leads to cleaner code overall. Please see the documentation of
`
nifty7.random
` for details.
Interface Change for from_random and OuterProduct
-------------------------------------------------
...
...
src/field.py
View file @
04c52785
...
...
@@ -136,6 +136,8 @@ class Field(Operator):
The domain of the output random Field.
dtype : type
The datatype of the output random Field.
If the datatype is complex, each real and imaginary part
have variance 1
Returns
-------
...
...
src/multi_field.py
View file @
04c52785
...
...
@@ -112,6 +112,8 @@ class MultiField(Operator):
The domain of the output random Field.
dtype : type
The datatype of the output random Field.
If the datatype is complex, each real an imaginary part have
variance 1.
Returns
-------
...
...
src/operators/energy_operators.py
View file @
04c52785
...
...
@@ -237,6 +237,13 @@ class GaussianEnergy(EnergyOperator):
domain : Domain, DomainTuple, tuple of Domain or MultiDomain
Operator domain. By default it is inferred from `mean` or
`covariance` if specified
sampling_dtype : type
Here one can specify whether the distribution is a complex Gaussian or
not. Note that for a complex Gaussian the inverse_covariance is
.. math ::
(<ff^dagger>)^{-1}_P(f)/2,
where the additional factor of 2 is necessary because the
domain of s has double as many dimensions as in the real case.
Note
----
...
...
src/random.py
View file @
04c52785
...
...
@@ -238,8 +238,8 @@ class Random(object):
raise
TypeError
(
"mean must not be complex for a real result field"
)
if
np
.
issubdtype
(
dtype
,
np
.
complexfloating
):
x
=
np
.
empty
(
shape
,
dtype
=
dtype
)
x
.
real
=
_rng
[
-
1
].
normal
(
mean
.
real
,
std
*
np
.
sqrt
(
0.5
)
,
shape
)
x
.
imag
=
_rng
[
-
1
].
normal
(
mean
.
imag
,
std
*
np
.
sqrt
(
0.5
)
,
shape
)
x
.
real
=
_rng
[
-
1
].
normal
(
mean
.
real
,
std
,
shape
)
x
.
imag
=
_rng
[
-
1
].
normal
(
mean
.
imag
,
std
,
shape
)
else
:
x
=
_rng
[
-
1
].
normal
(
mean
,
std
,
shape
).
astype
(
dtype
,
copy
=
False
)
return
x
...
...
src/sugar.py
View file @
04c52785
...
...
@@ -266,6 +266,8 @@ def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
The random distribution to use.
dtype : type
data type of the output field (e.g. numpy.float64)
If the datatype is complex, each real an imaginary part have
variance 1.
**kwargs : additional parameters for the random distribution
('mean' and 'std' for 'normal', 'low' and 'high' for 'uniform')
...
...
test/test_operators/test_complex_consistency.py
0 → 100644
View file @
04c52785
# 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-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import
numpy
as
np
import
nifty7
as
ift
Nsamp
=
20000
np
.
random
.
seed
(
42
)
def
_to_array
(
d
):
if
isinstance
(
d
,
np
.
ndarray
):
return
d
assert
isinstance
(
d
,
dict
)
return
np
.
concatenate
(
list
(
d
.
values
()))
def
test_GaussianEnergy
():
sp
=
ift
.
UnstructuredDomain
(
Nsamp
)
S
=
ift
.
ScalingOperator
(
sp
,
1.
)
samp
=
S
.
draw_sample_with_dtype
(
dtype
=
np
.
complex128
)
real_std
=
np
.
std
(
samp
.
val
.
real
)
imag_std
=
np
.
std
(
samp
.
val
.
imag
)
np
.
testing
.
assert_allclose
(
real_std
,
imag_std
,
atol
=
5.
/
np
.
sqrt
(
Nsamp
))
sp1
=
ift
.
UnstructuredDomain
(
1
)
mean
=
ift
.
full
(
sp1
,
0.
)
real_icov
=
ift
.
ScalingOperator
(
sp1
,
real_std
**
(
-
2
))
imag_icov
=
ift
.
ScalingOperator
(
sp1
,
imag_std
**
(
-
2
))
real_energy
=
ift
.
GaussianEnergy
(
mean
,
inverse_covariance
=
real_icov
)
imag_energy
=
ift
.
GaussianEnergy
(
mean
,
inverse_covariance
=
imag_icov
)
icov
=
ift
.
ScalingOperator
(
sp1
,
1.
)
complex_energy
=
ift
.
GaussianEnergy
(
mean
+
0.j
,
inverse_covariance
=
icov
)
for
i
in
range
(
min
(
10
,
Nsamp
)):
fld
=
ift
.
full
(
sp1
,
samp
.
val
[
i
])
val1
=
(
real_energy
(
fld
.
real
)
+
imag_energy
(
fld
.
imag
)).
val
val2
=
complex_energy
(
fld
).
val
np
.
testing
.
assert_allclose
(
val1
,
val2
,
atol
=
10.
/
np
.
sqrt
(
Nsamp
))
test/test_operators/test_fisher_metric.py
0 → 100644
View file @
04c52785
# 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-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import
numpy
as
np
import
pytest
import
nifty7
as
ift
from
..common
import
list2fixture
,
setup_function
,
teardown_function
spaces
=
[
ift
.
GLSpace
(
5
),
ift
.
MultiDomain
.
make
({
''
:
ift
.
RGSpace
(
5
,
distances
=
.
789
)}),
(
ift
.
RGSpace
(
3
,
distances
=
.
789
),
ift
.
UnstructuredDomain
(
2
))]
pmp
=
pytest
.
mark
.
parametrize
field
=
list2fixture
([
ift
.
from_random
(
sp
,
'normal'
)
for
sp
in
spaces
]
+
[
ift
.
from_random
(
sp
,
'normal'
,
dtype
=
np
.
complex128
)
for
sp
in
spaces
])
dtype
=
list2fixture
([
np
.
float64
,
np
.
complex128
])
Nsamp
=
2000
np
.
random
.
seed
(
42
)
def
_to_array
(
d
):
if
isinstance
(
d
,
np
.
ndarray
):
return
d
assert
isinstance
(
d
,
dict
)
return
np
.
concatenate
(
list
(
d
.
values
()))
def
_complex2real
(
sp
):
tup
=
tuple
([
d
for
d
in
sp
])
rsp
=
ift
.
DomainTuple
.
make
((
ift
.
UnstructuredDomain
(
2
),)
+
tup
)
rl
=
ift
.
DomainTupleFieldInserter
(
rsp
,
0
,
(
0
,))
im
=
ift
.
DomainTupleFieldInserter
(
rsp
,
0
,
(
1
,))
x
=
ift
.
ScalingOperator
(
sp
,
1
)
return
rl
(
x
.
real
)
+
im
(
x
.
imag
)
def
test_complex2real
():
sp
=
ift
.
UnstructuredDomain
(
3
)
op
=
_complex2real
(
ift
.
makeDomain
(
sp
))
f
=
ift
.
from_random
(
op
.
domain
,
'normal'
,
dtype
=
np
.
complex128
)
assert
np
.
all
((
f
==
op
.
adjoint_times
(
op
(
f
))).
val
)
assert
op
(
f
).
dtype
==
np
.
float64
f
=
ift
.
from_random
(
op
.
target
,
'normal'
)
assert
np
.
all
((
f
==
op
(
op
.
adjoint_times
(
f
))).
val
)
def
energy_tester
(
pos
,
get_noisy_data
,
energy_initializer
):
if
isinstance
(
pos
,
ift
.
Field
):
if
np
.
issubdtype
(
pos
.
dtype
,
np
.
complexfloating
):
op
=
_complex2real
(
pos
.
domain
)
else
:
op
=
ift
.
ScalingOperator
(
pos
.
domain
,
1.
)
else
:
ops
=
[]
for
k
,
dom
in
pos
.
domain
.
items
():
if
np
.
issubdtype
(
pos
[
k
].
dtype
,
np
.
complexfloating
):
ops
.
append
(
_complex2real
(
dom
).
ducktape
(
k
).
ducktape_left
(
k
))
else
:
FA
=
ift
.
FieldAdapter
(
dom
,
k
)
ops
.
append
(
FA
.
adjoint
@
FA
)
realizer
=
ift
.
utilities
.
my_sum
(
ops
)
from
nifty7.operator_spectrum
import
_DomRemover
flattener
=
_DomRemover
(
realizer
.
target
)
op
=
flattener
@
realizer
npos
=
op
(
pos
)
nget_noisy_data
=
lambda
mean
:
get_noisy_data
(
op
.
adjoint_times
(
mean
))
nenergy_initializer
=
lambda
mean
:
energy_initializer
(
mean
)
@
op
.
adjoint
_actual_energy_tester
(
npos
,
nget_noisy_data
,
nenergy_initializer
)
def
_actual_energy_tester
(
pos
,
get_noisy_data
,
energy_initializer
):
domain
=
pos
.
domain
test_vec
=
ift
.
from_random
(
domain
,
'normal'
)
results
=
[]
lin
=
ift
.
Linearization
.
make_var
(
pos
)
for
i
in
range
(
Nsamp
):
data
=
get_noisy_data
(
pos
)
energy
=
energy_initializer
(
data
)
grad
=
energy
(
lin
).
jac
.
adjoint
(
ift
.
full
(
energy
.
target
,
1.
))
results
.
append
(
_to_array
((
grad
*
grad
.
s_vdot
(
test_vec
)).
val
))
print
(
energy
)
print
(
grad
)
res
=
np
.
mean
(
np
.
array
(
results
),
axis
=
0
)
std
=
np
.
std
(
np
.
array
(
results
),
axis
=
0
)
/
np
.
sqrt
(
Nsamp
)
energy
=
energy_initializer
(
data
)
lin
=
ift
.
Linearization
.
make_var
(
pos
,
want_metric
=
True
)
res2
=
_to_array
(
energy
(
lin
).
metric
(
test_vec
).
val
)
np
.
testing
.
assert_allclose
(
res
/
std
,
res2
/
std
,
atol
=
6
)
def
test_GaussianEnergy
(
field
):
dtype
=
field
.
dtype
icov
=
ift
.
from_random
(
field
.
domain
,
'normal'
)
**
2
icov
=
ift
.
makeOp
(
icov
)
get_noisy_data
=
lambda
mean
:
mean
+
icov
.
draw_sample_with_dtype
(
from_inverse
=
True
,
dtype
=
dtype
)
E_init
=
lambda
data
:
ift
.
GaussianEnergy
(
mean
=
data
,
inverse_covariance
=
icov
)
energy_tester
(
field
,
get_noisy_data
,
E_init
)
def
test_PoissonEnergy
(
field
):
if
not
isinstance
(
field
,
ift
.
Field
):
pytest
.
skip
(
"MultiField Poisson energy not supported"
)
if
np
.
iscomplexobj
(
field
.
val
):
pytest
.
skip
(
"Poisson energy not defined for complex flux"
)
get_noisy_data
=
lambda
mean
:
ift
.
makeField
(
mean
.
domain
,
np
.
random
.
poisson
(
mean
.
val
))
# Make rate positive and high enough to avoid bad statistic
lam
=
10
*
(
field
**
2
).
clip
(
0.1
,
None
)
E_init
=
lambda
data
:
ift
.
PoissonianEnergy
(
data
)
energy_tester
(
lam
,
get_noisy_data
,
E_init
)
def
test_VariableCovarianceGaussianEnergy
(
dtype
):
dom
=
ift
.
UnstructuredDomain
(
3
)
res
=
ift
.
from_random
(
dom
,
'normal'
,
dtype
=
dtype
)
ivar
=
ift
.
from_random
(
dom
,
'normal'
)
**
2
+
4.
mf
=
ift
.
MultiField
.
from_dict
({
'res'
:
res
,
'ivar'
:
ivar
})
energy
=
ift
.
VariableCovarianceGaussianEnergy
(
dom
,
'res'
,
'ivar'
,
dtype
)
def
get_noisy_data
(
mean
):
samp
=
ift
.
from_random
(
dom
,
'normal'
,
dtype
)
samp
=
samp
/
mean
[
'ivar'
].
sqrt
()
return
samp
+
mean
[
'res'
]
def
E_init
(
data
):
adder
=
ift
.
Adder
(
ift
.
MultiField
.
from_dict
({
'res'
:
data
}),
neg
=
True
)
return
energy
.
partial_insert
(
adder
)
energy_tester
(
mf
,
get_noisy_data
,
E_init
)
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new 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