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
df7b30c9
Commit
df7b30c9
authored
Jul 14, 2018
by
Martin Reinecke
Browse files
merge
parents
bedbdb64
e8105d79
Changes
20
Hide whitespace changes
Inline
Side-by-side
demos/getting_started_1.py
View file @
df7b30c9
...
...
@@ -115,5 +115,5 @@ if __name__ == '__main__':
ift
.
plot
(
mask_to_nan
(
mask
,
(
GR
*
Mask
).
adjoint
(
data
)),
title
=
'Data'
)
ift
.
plot
(
HT
(
m
),
title
=
'Reconstruction'
)
ift
.
plot
(
mask_to_nan
(
mask
,
HT
(
m
-
MOCK_SIGNAL
)))
ift
.
plot_finish
(
nx
=
4
,
ny
=
1
,
xsize
=
2
0
,
ysize
=
4
,
ift
.
plot_finish
(
nx
=
2
,
ny
=
2
,
xsize
=
1
0
,
ysize
=
10
,
title
=
"getting_started_1"
)
nifty5/__init__.py
View file @
df7b30c9
...
...
@@ -32,7 +32,9 @@ from .operators.domain_distributor import DomainDistributor
from
.operators.endomorphic_operator
import
EndomorphicOperator
from
.operators.exp_transform
import
ExpTransform
from
.operators.fft_operator
import
FFTOperator
from
.operators.fft_smoothing_operator
import
FFTSmoothingOperator
from
.operators.field_zero_padder
import
FieldZeroPadder
from
.operators.hartley_operator
import
HartleyOperator
from
.operators.harmonic_smoothing_operator
import
HarmonicSmoothingOperator
from
.operators.geometry_remover
import
GeometryRemover
from
.operators.harmonic_transform_operator
import
HarmonicTransformOperator
from
.operators.inversion_enabler
import
InversionEnabler
...
...
nifty5/domain_tuple.py
View file @
df7b30c9
...
...
@@ -52,7 +52,7 @@ class DomainTuple(object):
nax
=
len
(
thing
.
shape
)
res
[
idx
]
=
tuple
(
range
(
i
,
i
+
nax
))
i
+=
nax
return
res
return
tuple
(
res
)
@
staticmethod
def
make
(
domain
):
...
...
nifty5/library/correlated_fields.py
View file @
df7b30c9
...
...
@@ -25,7 +25,7 @@ from ..models.local_nonlinearity import PointwiseExponential
from
..models.variable
import
Variable
from
..multi.multi_field
import
MultiField
from
..operators.domain_distributor
import
DomainDistributor
from
..operators.
fft
_operator
import
FFT
Operator
from
..operators.
hartley
_operator
import
Hartley
Operator
from
..operators.harmonic_transform_operator
import
HarmonicTransformOperator
from
..operators.power_distributor
import
PowerDistributor
...
...
@@ -41,9 +41,10 @@ def make_correlated_field(s_space, amplitude_model):
amplitude_model : model for correlation structure
'''
h_space
=
s_space
.
get_default_codomain
()
ht
=
FFT
Operator
(
h_space
,
s_space
)
ht
=
Hartley
Operator
(
h_space
,
s_space
)
p_space
=
amplitude_model
.
value
.
domain
[
0
]
power_distributor
=
PowerDistributor
(
h_space
,
p_space
)
# FIXME Remove tau and phi stuff from here. Should not be necessary
position
=
MultiField
.
from_dict
({
'xi'
:
Field
.
from_random
(
'normal'
,
h_space
),
'tau'
:
amplitude_model
.
position
[
'tau'
],
...
...
nifty5/models/constant.py
View file @
df7b30c9
...
...
@@ -37,9 +37,8 @@ class Constant(Model):
-----
Since there is no model-function associated:
- Position has no influence on value.
- The
re is no
Jacobian.
- The Jacobian
is a null matrix
.
"""
# TODO Remove position
def
__init__
(
self
,
position
,
constant
):
super
(
Constant
,
self
).
__init__
(
position
)
self
.
_constant
=
constant
...
...
nifty5/multi/multi_field.py
View file @
df7b30c9
...
...
@@ -242,7 +242,7 @@ for op in ["__sub__", "__rsub__",
if
self
.
_domain
is
not
other
.
_domain
:
raise
ValueError
(
"domain mismatch"
)
val
=
tuple
(
getattr
(
v1
,
op
)(
v2
)
for
v1
,
v2
in
zip
(
self
.
_val
,
other
.
_val
))
for
v1
,
v2
in
zip
(
self
.
_val
,
other
.
_val
))
else
:
val
=
tuple
(
getattr
(
v1
,
op
)(
other
)
for
v1
in
self
.
_val
)
return
MultiField
(
self
.
_domain
,
val
)
...
...
nifty5/operators/domain_distributor.py
View file @
df7b30c9
...
...
@@ -25,26 +25,16 @@ from ..compat import *
from
..domain_tuple
import
DomainTuple
from
..field
import
Field
from
.linear_operator
import
LinearOperator
from
..
import
utilities
# MR FIXME: this needs to be rewritten in a generic fashion
class
DomainDistributor
(
LinearOperator
):
def
__init__
(
self
,
target
,
axis
):
if
dobj
.
ntask
>
1
:
raise
NotImplementedError
(
'UpProj class does not support MPI.'
)
assert
len
(
target
)
==
2
assert
axis
in
[
0
,
1
]
if
axis
==
0
:
domain
=
target
[
1
]
self
.
_size
=
target
[
0
].
size
else
:
domain
=
target
[
0
]
self
.
_size
=
target
[
1
].
size
self
.
_axis
=
axis
self
.
_domain
=
DomainTuple
.
make
(
domain
)
def
__init__
(
self
,
target
,
spaces
):
self
.
_target
=
DomainTuple
.
make
(
target
)
self
.
_spaces
=
utilities
.
parse_spaces
(
spaces
,
len
(
self
.
_target
))
self
.
_domain
=
[
tgt
for
i
,
tgt
in
enumerate
(
self
.
_target
)
if
i
in
self
.
_spaces
]
self
.
_domain
=
DomainTuple
.
make
(
self
.
_domain
)
@
property
def
domain
(
self
):
...
...
@@ -57,23 +47,16 @@ class DomainDistributor(LinearOperator):
def
apply
(
self
,
x
,
mode
):
self
.
_check_input
(
x
,
mode
)
if
mode
==
self
.
TIMES
:
x
=
x
.
local_data
otherDirection
=
np
.
ones
(
self
.
_size
)
if
self
.
_axis
==
0
:
res
=
np
.
outer
(
otherDirection
,
x
)
else
:
res
=
np
.
outer
(
x
,
otherDirection
)
res
=
res
.
reshape
(
dobj
.
local_shape
(
self
.
target
.
shape
))
return
Field
.
from_local_data
(
self
.
target
,
res
)
ldat
=
x
.
local_data
if
0
in
self
.
_spaces
else
x
.
to_global_data
()
shp
=
[]
for
i
,
tgt
in
enumerate
(
self
.
_target
):
tmp
=
tgt
.
shape
if
i
>
0
else
tgt
.
local_shape
shp
+=
tmp
if
i
in
self
.
_spaces
else
(
1
,)
*
len
(
tgt
.
shape
)
ldat
=
np
.
broadcast_to
(
ldat
.
reshape
(
shp
),
self
.
_target
.
local_shape
)
return
Field
.
from_local_data
(
self
.
_target
,
ldat
)
else
:
if
self
.
_axis
==
0
:
x
=
x
.
local_data
.
reshape
(
self
.
_size
,
-
1
)
res
=
np
.
sum
(
x
,
axis
=
0
)
else
:
x
=
x
.
local_data
.
reshape
(
-
1
,
self
.
_size
)
res
=
np
.
sum
(
x
,
axis
=
1
)
res
=
res
.
reshape
(
dobj
.
local_shape
(
self
.
domain
.
shape
))
return
Field
.
from_local_data
(
self
.
domain
,
res
)
return
x
.
sum
([
s
for
s
in
range
(
len
(
x
.
domain
))
if
s
not
in
self
.
_spaces
])
@
property
def
capability
(
self
):
...
...
nifty5/operators/exp_transform.py
View file @
df7b30c9
...
...
@@ -27,20 +27,23 @@ from ..domains.power_space import PowerSpace
from
..domains.rg_space
import
RGSpace
from
..field
import
Field
from
.linear_operator
import
LinearOperator
from
..
import
utilities
class
ExpTransform
(
LinearOperator
):
def
__init__
(
self
,
target
,
dof
):
if
not
((
isinstance
(
target
,
RGSpace
)
and
target
.
harmonic
)
or
isinstance
(
target
,
PowerSpace
)):
def
__init__
(
self
,
target
,
dof
,
space
=
0
):
self
.
_target
=
DomainTuple
.
make
(
target
)
self
.
_space
=
utilities
.
infer_space
(
self
.
_target
,
space
)
tgt
=
self
.
_target
[
self
.
_space
]
if
not
((
isinstance
(
tgt
,
RGSpace
)
and
tgt
.
harmonic
)
or
isinstance
(
tgt
,
PowerSpace
)):
raise
ValueError
(
"Target must be a harmonic RGSpace or a power space."
)
if
np
.
isscalar
(
dof
):
dof
=
np
.
full
(
len
(
t
arge
t
.
shape
),
int
(
dof
),
dtype
=
np
.
int
)
dof
=
np
.
full
(
len
(
t
g
t
.
shape
),
int
(
dof
),
dtype
=
np
.
int
)
dof
=
np
.
array
(
dof
)
ndim
=
len
(
t
arge
t
.
shape
)
ndim
=
len
(
t
g
t
.
shape
)
t_mins
=
np
.
empty
(
ndim
)
bindistances
=
np
.
empty
(
ndim
)
...
...
@@ -48,12 +51,12 @@ class ExpTransform(LinearOperator):
self
.
_frac
=
[
None
]
*
ndim
for
i
in
range
(
ndim
):
if
isinstance
(
t
arge
t
,
RGSpace
):
rng
=
np
.
arange
(
t
arge
t
.
shape
[
i
])
tmp
=
np
.
minimum
(
rng
,
t
arge
t
.
shape
[
i
]
+
1
-
rng
)
k_array
=
tmp
*
t
arge
t
.
distances
[
i
]
if
isinstance
(
t
g
t
,
RGSpace
):
rng
=
np
.
arange
(
t
g
t
.
shape
[
i
])
tmp
=
np
.
minimum
(
rng
,
t
g
t
.
shape
[
i
]
+
1
-
rng
)
k_array
=
tmp
*
t
g
t
.
distances
[
i
]
else
:
k_array
=
t
arge
t
.
k_lengths
k_array
=
t
g
t
.
k_lengths
# avoid taking log of first entry
log_k_array
=
np
.
log
(
k_array
[
1
:])
...
...
@@ -77,8 +80,9 @@ class ExpTransform(LinearOperator):
from
..domains.log_rg_space
import
LogRGSpace
log_space
=
LogRGSpace
(
2
*
dof
+
1
,
bindistances
,
t_mins
,
harmonic
=
False
)
self
.
_target
=
DomainTuple
.
make
(
target
)
self
.
_domain
=
DomainTuple
.
make
(
log_space
)
self
.
_domain
=
[
dom
for
dom
in
self
.
_target
]
self
.
_domain
[
self
.
_space
]
=
log_space
self
.
_domain
=
DomainTuple
.
make
(
self
.
_domain
)
@
property
def
domain
(
self
):
...
...
@@ -94,9 +98,10 @@ class ExpTransform(LinearOperator):
ax
=
dobj
.
distaxis
(
x
)
ndim
=
len
(
self
.
target
.
shape
)
curshp
=
list
(
self
.
_dom
(
mode
).
shape
)
for
d
in
range
(
ndim
):
d0
=
self
.
_target
.
axes
[
self
.
_space
][
0
]
for
d
in
self
.
_target
.
axes
[
self
.
_space
]:
idx
=
(
slice
(
None
,),)
*
d
wgt
=
self
.
_frac
[
d
].
reshape
((
1
,)
*
d
+
(
-
1
,)
+
(
1
,)
*
(
ndim
-
d
-
1
))
wgt
=
self
.
_frac
[
d
-
d0
].
reshape
((
1
,)
*
d
+
(
-
1
,)
+
(
1
,)
*
(
ndim
-
d
-
1
))
if
d
==
ax
:
x
=
dobj
.
redistribute
(
x
,
nodist
=
(
ax
,))
...
...
@@ -107,11 +112,11 @@ class ExpTransform(LinearOperator):
shp
=
list
(
x
.
shape
)
shp
[
d
]
=
self
.
_tgt
(
mode
).
shape
[
d
]
xnew
=
np
.
zeros
(
shp
,
dtype
=
x
.
dtype
)
np
.
add
.
at
(
xnew
,
idx
+
(
self
.
_bindex
[
d
],),
x
*
(
1.
-
wgt
))
np
.
add
.
at
(
xnew
,
idx
+
(
self
.
_bindex
[
d
]
+
1
,),
x
*
wgt
)
np
.
add
.
at
(
xnew
,
idx
+
(
self
.
_bindex
[
d
-
d0
],),
x
*
(
1.
-
wgt
))
np
.
add
.
at
(
xnew
,
idx
+
(
self
.
_bindex
[
d
-
d0
]
+
1
,),
x
*
wgt
)
else
:
# TIMES
xnew
=
x
[
idx
+
(
self
.
_bindex
[
d
],)]
*
(
1.
-
wgt
)
xnew
+=
x
[
idx
+
(
self
.
_bindex
[
d
]
+
1
,)]
*
wgt
xnew
=
x
[
idx
+
(
self
.
_bindex
[
d
-
d0
],)]
*
(
1.
-
wgt
)
xnew
+=
x
[
idx
+
(
self
.
_bindex
[
d
-
d0
]
+
1
,)]
*
wgt
curshp
[
d
]
=
self
.
_tgt
(
mode
).
shape
[
d
]
x
=
dobj
.
from_local_data
(
curshp
,
xnew
,
distaxis
=
curax
)
...
...
nifty5/operators/fft_operator.py
View file @
df7b30c9
...
...
@@ -43,6 +43,13 @@ class FFTOperator(LinearOperator):
The index of the subdomain on which the operator should act
If None, it is set to 0 if `domain` contains exactly one space.
`domain[space]` must be an RGSpace.
Notes
-----
This operator performs full FFTs, which implies that its output field will
always have complex type, regardless of the type of the input field.
If a real field is desired after a forward/backward transform couple, it
must be manually cast to real.
"""
def
__init__
(
self
,
domain
,
target
=
None
,
space
=
None
):
...
...
@@ -67,41 +74,35 @@ class FFTOperator(LinearOperator):
utilities
.
fft_prep
()
def
apply
(
self
,
x
,
mode
):
from
pyfftw.interfaces.numpy_fft
import
fftn
,
ifftn
self
.
_check_input
(
x
,
mode
)
if
np
.
issubdtype
(
x
.
dtype
,
np
.
complexfloating
):
return
(
self
.
_apply_cartesian
(
x
.
real
,
mode
)
+
1j
*
self
.
_apply_cartesian
(
x
.
imag
,
mode
))
ncells
=
x
.
domain
[
self
.
_space
].
size
if
x
.
domain
[
self
.
_space
].
harmonic
:
# harmonic -> position
func
=
fftn
fct
=
1.
else
:
return
self
.
_apply_cartesian
(
x
,
mode
)
def
_apply_cartesian
(
self
,
x
,
mode
):
func
=
ifftn
fct
=
ncells
axes
=
x
.
domain
.
axes
[
self
.
_space
]
tdom
=
self
.
_tgt
(
mode
)
oldax
=
dobj
.
distaxis
(
x
.
val
)
if
oldax
not
in
axes
:
# straightforward, no redistribution needed
ldat
=
x
.
local_data
ldat
=
utilities
.
hartley
(
ldat
,
axes
=
axes
)
ldat
=
func
(
ldat
,
axes
=
axes
)
tmp
=
dobj
.
from_local_data
(
x
.
val
.
shape
,
ldat
,
distaxis
=
oldax
)
elif
len
(
axes
)
<
len
(
x
.
shape
)
or
len
(
axes
)
==
1
:
# we can use one
Hartley
pass in between the redistributions
# we can use one
FFT
pass in between the redistributions
tmp
=
dobj
.
redistribute
(
x
.
val
,
nodist
=
axes
)
newax
=
dobj
.
distaxis
(
tmp
)
ldat
=
dobj
.
local_data
(
tmp
)
ldat
=
utilities
.
hartley
(
ldat
,
axes
=
axes
)
ldat
=
func
(
ldat
,
axes
=
axes
)
tmp
=
dobj
.
from_local_data
(
tmp
.
shape
,
ldat
,
distaxis
=
newax
)
tmp
=
dobj
.
redistribute
(
tmp
,
dist
=
oldax
)
else
:
# two separate, full FFTs needed
# ideal strategy for the moment would be:
# - do real-to-complex FFT on all local axes
# - fill up array
# - redistribute array
# - do complex-to-complex FFT on remaining axis
# - add re+im
# - redistribute back
else
:
# two separate FFTs needed
rem_axes
=
tuple
(
i
for
i
in
axes
if
i
!=
oldax
)
tmp
=
x
.
val
ldat
=
dobj
.
local_data
(
tmp
)
ldat
=
utilities
.
my_fftn_r2
c
(
ldat
,
axes
=
rem_axes
)
ldat
=
fun
c
(
ldat
,
axes
=
rem_axes
)
if
oldax
!=
0
:
raise
ValueError
(
"bad distribution"
)
ldat2
=
ldat
.
reshape
((
ldat
.
shape
[
0
],
...
...
@@ -110,17 +111,16 @@ class FFTOperator(LinearOperator):
tmp
=
dobj
.
from_local_data
(
shp2d
,
ldat2
,
distaxis
=
0
)
tmp
=
dobj
.
transpose
(
tmp
)
ldat2
=
dobj
.
local_data
(
tmp
)
ldat2
=
utilities
.
my_fftn
(
ldat2
,
axes
=
(
1
,))
ldat2
=
ldat2
.
real
+
ldat2
.
imag
ldat2
=
func
(
ldat2
,
axes
=
(
1
,))
tmp
=
dobj
.
from_local_data
(
tmp
.
shape
,
ldat2
,
distaxis
=
0
)
tmp
=
dobj
.
transpose
(
tmp
)
ldat2
=
dobj
.
local_data
(
tmp
).
reshape
(
ldat
.
shape
)
tmp
=
dobj
.
from_local_data
(
x
.
val
.
shape
,
ldat2
,
distaxis
=
0
)
Tval
=
Field
(
tdom
,
tmp
)
if
mode
&
(
LinearOperator
.
TIMES
|
LinearOperator
.
ADJOINT_TIMES
):
fct
=
self
.
_domain
[
self
.
_space
].
scalar_dvol
fct
*
=
self
.
_domain
[
self
.
_space
].
scalar_dvol
else
:
fct
=
self
.
_target
[
self
.
_space
].
scalar_dvol
fct
*
=
self
.
_target
[
self
.
_space
].
scalar_dvol
return
Tval
if
fct
==
1
else
Tval
*
fct
@
property
...
...
nifty5/operators/field_zero_padder.py
0 → 100644
View file @
df7b30c9
from
__future__
import
absolute_import
,
division
,
print_function
import
numpy
as
np
from
..
import
dobj
from
..compat
import
*
from
..domain_tuple
import
DomainTuple
from
..domains.rg_space
import
RGSpace
from
..field
import
Field
from
.linear_operator
import
LinearOperator
from
..
import
utilities
class
FieldZeroPadder
(
LinearOperator
):
def
__init__
(
self
,
domain
,
factor
,
space
=
0
):
super
(
FieldZeroPadder
,
self
).
__init__
()
self
.
_domain
=
DomainTuple
.
make
(
domain
)
self
.
_space
=
utilities
.
infer_space
(
self
.
_domain
,
space
)
dom
=
self
.
_domain
[
self
.
_space
]
if
not
isinstance
(
dom
,
RGSpace
):
raise
TypeError
(
"RGSpace required"
)
if
not
len
(
dom
.
shape
)
==
1
:
raise
TypeError
(
"RGSpace must be one-dimensional"
)
if
dom
.
harmonic
:
raise
TypeError
(
"RGSpace must not be harmonic"
)
tgt
=
RGSpace
((
int
(
factor
*
dom
.
shape
[
0
]),),
dom
.
distances
)
self
.
_target
=
list
(
self
.
_domain
)
self
.
_target
[
self
.
_space
]
=
tgt
self
.
_target
=
DomainTuple
.
make
(
self
.
_target
)
@
property
def
domain
(
self
):
return
self
.
_domain
@
property
def
target
(
self
):
return
self
.
_target
@
property
def
capability
(
self
):
return
self
.
TIMES
|
self
.
ADJOINT_TIMES
def
apply
(
self
,
x
,
mode
):
self
.
_check_input
(
x
,
mode
)
x
=
x
.
val
dax
=
dobj
.
distaxis
(
x
)
shp_in
=
x
.
shape
shp_out
=
self
.
_tgt
(
mode
).
shape
ax
=
self
.
_target
.
axes
[
self
.
_space
][
0
]
if
dax
==
ax
:
x
=
dobj
.
redistribute
(
x
,
nodist
=
(
ax
,))
curax
=
dobj
.
distaxis
(
x
)
if
mode
==
self
.
ADJOINT_TIMES
:
newarr
=
np
.
empty
(
dobj
.
local_shape
(
shp_out
,
curax
),
dtype
=
x
.
dtype
)
newarr
[()]
=
dobj
.
local_data
(
x
)[(
slice
(
None
),)
*
ax
+
(
slice
(
0
,
shp_out
[
ax
]),)]
else
:
newarr
=
np
.
zeros
(
dobj
.
local_shape
(
shp_out
,
curax
),
dtype
=
x
.
dtype
)
newarr
[(
slice
(
None
),)
*
ax
+
(
slice
(
0
,
shp_in
[
ax
]),)]
=
dobj
.
local_data
(
x
)
newarr
=
dobj
.
from_local_data
(
shp_out
,
newarr
,
distaxis
=
curax
)
if
dax
==
ax
:
newarr
=
dobj
.
redistribute
(
newarr
,
dist
=
ax
)
return
Field
(
self
.
_tgt
(
mode
),
val
=
newarr
)
nifty5/operators/
fft
_smoothing_operator.py
→
nifty5/operators/
harmonic
_smoothing_operator.py
View file @
df7b30c9
...
...
@@ -22,11 +22,11 @@ from ..compat import *
from
..domain_tuple
import
DomainTuple
from
..utilities
import
infer_space
from
.diagonal_operator
import
DiagonalOperator
from
.
fft
_operator
import
FFT
Operator
from
.
hartley
_operator
import
Hartley
Operator
from
.scaling_operator
import
ScalingOperator
def
FFT
SmoothingOperator
(
domain
,
sigma
,
space
=
None
):
def
Harmonic
SmoothingOperator
(
domain
,
sigma
,
space
=
None
):
""" This function returns an operator that carries out a smoothing with
a Gaussian kernel of width `sigma` on the part of `domain` given by
`space`.
...
...
@@ -59,12 +59,12 @@ def FFTSmoothingOperator(domain, sigma, space=None):
space
=
infer_space
(
domain
,
space
)
if
domain
[
space
].
harmonic
:
raise
TypeError
(
"domain must not be harmonic"
)
FFT
=
FFT
Operator
(
domain
,
space
=
space
)
codomain
=
FFT
.
domain
[
space
].
get_default_codomain
()
Hartley
=
Hartley
Operator
(
domain
,
space
=
space
)
codomain
=
Hartley
.
domain
[
space
].
get_default_codomain
()
kernel
=
codomain
.
get_k_length_array
()
smoother
=
codomain
.
get_fft_smoothing_kernel_function
(
sigma
)
kernel
=
smoother
(
kernel
)
ddom
=
list
(
domain
)
ddom
[
space
]
=
codomain
diag
=
DiagonalOperator
(
kernel
,
ddom
,
space
)
return
FFT
.
inverse
*
diag
*
FFT
return
Hartley
.
inverse
*
diag
*
Hartley
nifty5/operators/harmonic_transform_operator.py
View file @
df7b30c9
...
...
@@ -22,7 +22,7 @@ from .. import utilities
from
..compat
import
*
from
..domain_tuple
import
DomainTuple
from
..domains.rg_space
import
RGSpace
from
.
fft
_operator
import
FFT
Operator
from
.
hartley
_operator
import
Hartley
Operator
from
.linear_operator
import
LinearOperator
from
.sht_operator
import
SHTOperator
...
...
@@ -65,7 +65,7 @@ class HarmonicTransformOperator(LinearOperator):
raise
TypeError
(
"HarmonicTransformOperator only works on a harmonic space"
)
if
isinstance
(
hspc
,
RGSpace
):
self
.
_op
=
FFT
Operator
(
domain
,
target
,
space
)
self
.
_op
=
Hartley
Operator
(
domain
,
target
,
space
)
else
:
self
.
_op
=
SHTOperator
(
domain
,
target
,
space
)
...
...
nifty5/operators/hartley_operator.py
0 → 100644
View file @
df7b30c9
# 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-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from
__future__
import
absolute_import
,
division
,
print_function
import
numpy
as
np
from
..
import
dobj
,
utilities
from
..compat
import
*
from
..domain_tuple
import
DomainTuple
from
..domains.rg_space
import
RGSpace
from
..field
import
Field
from
.linear_operator
import
LinearOperator
class
HartleyOperator
(
LinearOperator
):
"""Transforms between a pair of position and harmonic RGSpaces.
Parameters
----------
domain: Domain, tuple of Domain or DomainTuple
The domain of the data that is input by "times" and output by
"adjoint_times".
target: Domain, optional
The target (sub-)domain of the transform operation.
If omitted, a domain will be chosen automatically.
space: int, optional
The index of the subdomain on which the operator should act
If None, it is set to 0 if `domain` contains exactly one space.
`domain[space]` must be an RGSpace.
Notes
-----
This operator always produces output fields with the same data type as
its input. This is achieved by performing so-called Hartley transforms
(https://en.wikipedia.org/wiki/Discrete_Hartley_transform).
For complex input fields, the operator will transform the real and
imaginary parts separately and use the results as real and imaginary parts
of the result field, respectivey.
In many contexts the Hartley transform is a perfect substitute for the
Fourier transform, but in some situations (e.g. convolution with a general,
non-symmetrc kernel, the full FFT must be used instead.
"""
def
__init__
(
self
,
domain
,
target
=
None
,
space
=
None
):
super
(
HartleyOperator
,
self
).
__init__
()
# Initialize domain and target
self
.
_domain
=
DomainTuple
.
make
(
domain
)
self
.
_space
=
utilities
.
infer_space
(
self
.
_domain
,
space
)
adom
=
self
.
_domain
[
self
.
_space
]
if
not
isinstance
(
adom
,
RGSpace
):
raise
TypeError
(
"HartleyOperator only works on RGSpaces"
)
if
target
is
None
:
target
=
adom
.
get_default_codomain
()
self
.
_target
=
[
dom
for
dom
in
self
.
_domain
]
self
.
_target
[
self
.
_space
]
=
target
self
.
_target
=
DomainTuple
.
make
(
self
.
_target
)
adom
.
check_codomain
(
target
)
target
.
check_codomain
(
adom
)
utilities
.
fft_prep
()
def
apply
(
self
,
x
,
mode
):
self
.
_check_input
(
x
,
mode
)
if
np
.
issubdtype
(
x
.
dtype
,
np
.
complexfloating
):
return
(
self
.
_apply_cartesian
(
x
.
real
,
mode
)
+
1j
*
self
.
_apply_cartesian
(
x
.
imag
,
mode
))
else
:
return
self
.
_apply_cartesian
(
x
,
mode
)
def
_apply_cartesian
(
self
,
x
,
mode
):
axes
=
x
.
domain
.
axes
[
self
.
_space
]
tdom
=
self
.
_tgt
(
mode
)
oldax
=
dobj
.
distaxis
(
x
.
val
)
if
oldax
not
in
axes
:
# straightforward, no redistribution needed
ldat
=
x
.
local_data
ldat
=
utilities
.
hartley
(
ldat
,
axes
=
axes
)
tmp
=
dobj
.
from_local_data
(
x
.
val
.
shape
,
ldat
,
distaxis
=
oldax
)
elif
len
(
axes
)
<
len
(
x
.
shape
)
or
len
(
axes
)
==
1
:
# we can use one Hartley pass in between the redistributions
tmp
=
dobj
.
redistribute
(
x
.
val
,
nodist
=
axes
)
newax
=
dobj
.
distaxis
(
tmp
)
ldat
=
dobj
.
local_data
(
tmp
)
ldat
=
utilities
.
hartley
(
ldat
,
axes
=
axes
)
tmp
=
dobj
.
from_local_data
(
tmp
.
shape
,
ldat
,
distaxis
=
newax
)
tmp
=
dobj
.
redistribute
(
tmp
,
dist
=
oldax
)
else
:
# two separate, full FFTs needed
# ideal strategy for the moment would be:
# - do real-to-complex FFT on all local axes
# - fill up array
# - redistribute array
# - do complex-to-complex FFT on remaining axis
# - add re+im
# - redistribute back
rem_axes
=
tuple
(
i
for
i
in
axes
if
i
!=
oldax
)
tmp
=
x
.
val
ldat
=
dobj
.
local_data
(
tmp
)
ldat
=
utilities
.
my_fftn_r2c
(
ldat
,
axes
=
rem_axes
)
if
oldax
!=
0
:
raise
ValueError
(
"bad distribution"
)
ldat2
=
ldat
.
reshape
((
ldat
.
shape
[
0
],
np
.
prod
(
ldat
.
shape
[
1
:])))
shp2d
=
(
x
.
val
.
shape
[
0
],
np
.
prod
(
x
.
val
.
shape
[
1
:]))
tmp
=
dobj
.
from_local_data
(
shp2d
,
ldat2
,
distaxis
=
0
)
tmp
=
dobj
.
transpose
(
tmp
)
ldat2
=
dobj
.
local_data
(
tmp
)
ldat2
=
utilities
.
my_fftn
(
ldat2
,
axes
=
(
1
,))
ldat2
=
ldat2
.
real
+
ldat2
.
imag
tmp
=
dobj
.
from_local_data
(
tmp
.
shape
,
ldat2
,
distaxis
=
0
)
tmp
=
dobj
.
transpose
(
tmp
)
ldat2
=
dobj
.
local_data
(
tmp
).
reshape
(
ldat
.
shape
)
tmp
=
dobj
.
from_local_data
(
x
.
val
.
shape
,
ldat2
,
distaxis
=
0
)
Tval
=
Field
(
tdom
,
tmp
)
if
mode
&
(
LinearOperator
.
TIMES
|
LinearOperator
.
ADJOINT_TIMES
):
fct
=
self
.
_domain
[
self
.
_space
].
scalar_dvol
else
:
fct
=
self
.
_target
[
self
.
_space
].
scalar_dvol
return
Tval
if
fct
==
1
else
Tval
*
fct
@
property
def
domain
(
self
):
return
self
.
_domain
@
property
def
target
(
self
):
return
self
.
_target
@
property
def
capability
(
self
):
return
self
.
_all_ops
nifty5/operators/linear_operator.py
View file @
df7b30c9
...
...
@@ -88,6 +88,7 @@ class LinearOperator(NiftyMetaBase()):
@
abc
.
abstractproperty