Commit b2f3e42a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge byebye_volume_factors branch

parents ca0d35fa 198eb267
Pipeline #23935 passed with stage
in 4 minutes and 46 seconds
...@@ -4,8 +4,8 @@ import numericalunits as nu ...@@ -4,8 +4,8 @@ import numericalunits as nu
if __name__ == "__main__": if __name__ == "__main__":
# In MPI mode, the random seed for numericalunits must be set by hand # In MPI mode, the random seed for numericalunits must be set by hand
nu.reset_units(43) #nu.reset_units(43)
dimensionality = 2 dimensionality = 1
np.random.seed(43) np.random.seed(43)
# Setting up variable parameters # Setting up variable parameters
...@@ -16,12 +16,13 @@ if __name__ == "__main__": ...@@ -16,12 +16,13 @@ if __name__ == "__main__":
field_sigma = 2. * nu.K field_sigma = 2. * nu.K
# smoothing length of response # smoothing length of response
response_sigma = 0.01*nu.m response_sigma = 0.01*nu.m
# The signal to noise ratio # The signal to noise ratio ***CURRENTLY BROKEN***
signal_to_noise = 0.7 signal_to_noise = 70.7
# note that field_variance**2 = a*k_0/4. for this analytic form of power # note that field_variance**2 = a*k_0/4. for this analytic form of power
# spectrum # spectrum
def power_spectrum(k): def power_spectrum(k):
#RL FIXME: signal_amplitude is not how much signal varies
cldim = correlation_length**(2*dimensionality) cldim = correlation_length**(2*dimensionality)
a = 4/(2*np.pi) * cldim * field_sigma**2 a = 4/(2*np.pi) * cldim * field_sigma**2
# to be integrated over spherical shells later on # to be integrated over spherical shells later on
...@@ -47,27 +48,34 @@ if __name__ == "__main__": ...@@ -47,27 +48,34 @@ if __name__ == "__main__":
mock_power = ift.PS_field(power_space, power_spectrum) mock_power = ift.PS_field(power_space, power_spectrum)
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True) mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
print mock_harmonic.val[0]/nu.K/(nu.m**dimensionality)
mock_signal = fft(mock_harmonic) mock_signal = fft(mock_harmonic)
print "msig",mock_signal.val[0:10]/nu.K
exposure = 1. sensitivity = (1./nu.m)**dimensionality/nu.K
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), R = ift.ResponseOperator(signal_space, sigma=(0.*response_sigma,),
exposure=(exposure,)) sensitivity=(sensitivity,))
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = R*fft R_harmonic = R*fft
noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality)
print noise_amplitude
N = ift.DiagonalOperator( N = ift.DiagonalOperator(
ift.Field.full(data_domain, ift.Field.full(data_domain, noise_amplitude**2))
mock_signal.var()/signal_to_noise).weight(1))
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=noise_amplitude, mean=0)
data = R(mock_signal) + noise data = R(mock_signal)
print data.val[5:10]
# Wiener filter data += noise
print data.val[5:10]
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R_harmonic.adjoint_times(N.inverse_times(data))
print "xx",j.val[0]*nu.K*(nu.m**dimensionality)
exit()
ctrl = ift.GradientNormController( ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-4/nu.K/(nu.m**(0.5*dimensionality))) verbose=True, tol_abs_gradnorm=1e-40/(nu.K*(nu.m**dimensionality)))
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter) S=S, N=N, R=R_harmonic, inverter=inverter)
...@@ -78,7 +86,7 @@ if __name__ == "__main__": ...@@ -78,7 +86,7 @@ if __name__ == "__main__":
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m) sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plot(ift.Field(sspace2, mock_signal.val)/nu.K, name="mock_signal.png") ift.plot(ift.Field(sspace2, mock_signal.val)/nu.K, name="mock_signal.png")
data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)/nu.K data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)
data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))/nu.K data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))
ift.plot(ift.Field(sspace2, val=data), name="data.png") ift.plot(ift.Field(sspace2, val=data), name="data.png")
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png") ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
from .scaling_operator import ScalingOperator
from .fft_operator import FFTOperator
from ..utilities import infer_space
from .diagonal_operator import DiagonalOperator
from .. import DomainTuple
def FFTSmoothingOperator(domain, sigma, space=None):
sigma = float(sigma)
if sigma < 0.:
raise ValueError("sigma must be nonnegative")
if sigma == 0.:
return ScalingOperator(1., domain)
domain = DomainTuple.make(domain)
space = infer_space(domain, space)
if domain[space].harmonic:
raise TypeError("domain must not be harmonic")
FFT = FFTOperator(domain, space=space)
codomain = FFT.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
...@@ -342,19 +342,10 @@ class Field(object): ...@@ -342,19 +342,10 @@ class Field(object):
spaces = utilities.parse_spaces(spaces, ndom) spaces = utilities.parse_spaces(spaces, ndom)
if len(spaces) == ndom: if len(spaces) == ndom:
tmp = self.scalar_weight(spaces) return dobj.vdot(self.val, x.val)
if tmp is None:
fct = 1.
y = self.weight(power=1)
else:
y = self
fct = tmp
return fct*dobj.vdot(y.val, x.val)
# If we arrive here, we have to do a partial dot product. # If we arrive here, we have to do a partial dot product.
# For the moment, do this the explicit, non-optimized way # For the moment, do this the explicit, non-optimized way
return (self.conjugate()*x).integrate(spaces=spaces) return (self.conjugate()*x).sum(spaces=spaces)
def norm(self): def norm(self):
""" Computes the L2-norm of the field values. """ Computes the L2-norm of the field values.
......
...@@ -74,10 +74,10 @@ class WienerFilterCurvature(EndomorphicOperator): ...@@ -74,10 +74,10 @@ class WienerFilterCurvature(EndomorphicOperator):
covariance. covariance.
""" """
power = sqrt(power_analyze(self.S.diagonal())) power = power_analyze(sqrt(self.S.diagonal()))
mock_signal = power_synthesize(power, real_signal=True) mock_signal = power_synthesize(power, real_signal=True)
noise = self.N.diagonal().weight(-1) noise = self.N.diagonal()
mock_noise = Field.from_random(random_type="normal", mock_noise = Field.from_random(random_type="normal",
domain=self.N.domain, dtype=noise.dtype) domain=self.N.domain, dtype=noise.dtype)
......
...@@ -82,7 +82,7 @@ class DOFProjectionOperator(LinearOperator): ...@@ -82,7 +82,7 @@ class DOFProjectionOperator(LinearOperator):
self._pshape = (presize, self._dofdex.size, postsize) self._pshape = (presize, self._dofdex.size, postsize)
def _times(self, x): def _times(self, x):
arr = dobj.local_data(x.weight(1).val) arr = dobj.local_data(x.val)
arr = arr.reshape(self._pshape) arr = arr.reshape(self._pshape)
oarr = np.zeros(self._hshape, dtype=x.dtype) oarr = np.zeros(self._hshape, dtype=x.dtype)
np.add.at(oarr, (slice(None), self._dofdex, slice(None)), arr) np.add.at(oarr, (slice(None), self._dofdex, slice(None)), arr)
...@@ -95,7 +95,7 @@ class DOFProjectionOperator(LinearOperator): ...@@ -95,7 +95,7 @@ class DOFProjectionOperator(LinearOperator):
res = Field(self._target, res = Field(self._target,
dobj.from_local_data(self._target.shape, oarr, dobj.from_local_data(self._target.shape, oarr,
dobj.default_distaxis())) dobj.default_distaxis()))
return res.weight(-1, spaces=self._space) return res
def _adjoint_times(self, x): def _adjoint_times(self, x):
res = Field.empty(self._domain, dtype=x.dtype) res = Field.empty(self._domain, dtype=x.dtype)
......
...@@ -19,9 +19,11 @@ ...@@ -19,9 +19,11 @@
import numpy as np import numpy as np
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..spaces.rg_space import RGSpace from ..spaces.rg_space import RGSpace
from ..utilities import infer_space
from .linear_operator import LinearOperator from .linear_operator import LinearOperator
from .fft_operator_support import RGRGTransformation, SphericalTransformation from .. import dobj
from .. import utilities
from ..field import Field
from ..spaces.gl_space import GLSpace
class FFTOperator(LinearOperator): class FFTOperator(LinearOperator):
...@@ -35,10 +37,11 @@ class FFTOperator(LinearOperator): ...@@ -35,10 +37,11 @@ class FFTOperator(LinearOperator):
Within a domain pair, both orderings are possible. Within a domain pair, both orderings are possible.
The operator provides a "times" and an "adjoint_times" operation. For RGSpaces, the operator provides the full set of operations.
For a pair of RGSpaces, the "adjoint_times" operation is equivalent to For the sphere-related domains, it only supports the transform from
"inverse_times"; for the sphere-related domains this is not the case, since harmonic to position space and its adjoint; if the operator domain is
the operator matrix is not square. harmonic, this will be times() and adjoint_times(), otherwise
inverse_times() and adjoint_inverse_times()
Parameters Parameters
---------- ----------
...@@ -63,35 +66,151 @@ class FFTOperator(LinearOperator): ...@@ -63,35 +66,151 @@ class FFTOperator(LinearOperator):
# Initialize domain and target # Initialize domain and target
self._domain = DomainTuple.make(domain) self._domain = DomainTuple.make(domain)
self._space = infer_space(self._domain, space) self._space = utilities.infer_space(self._domain, space)
adom = self.domain[self._space] adom = self._domain[self._space]
if target is None: if target is None:
target = adom.get_default_codomain() target = adom.get_default_codomain()
self._target = [dom for dom in self.domain] self._target = [dom for dom in self._domain]
self._target[self._space] = target self._target[self._space] = target
self._target = DomainTuple.make(self._target) self._target = DomainTuple.make(self._target)
adom.check_codomain(target) adom.check_codomain(target)
target.check_codomain(adom) target.check_codomain(adom)
if self._target[self._space].harmonic: if isinstance(adom, RGSpace):
pdom, hdom = (self._domain, self._target) self._applyfunc = self._apply_cartesian
self._capability = self._all_ops
import pyfftw
pyfftw.interfaces.cache.enable()
else: else:
pdom, hdom = (self._target, self._domain) from pyHealpix import sharpjob_d
if isinstance(pdom[self._space], RGSpace): self._applyfunc = self._apply_spherical
self._trafo = RGRGTransformation(pdom, hdom, self._space) hspc = adom if adom.harmonic else target
else: pspc = target if adom.harmonic else adom
self._trafo = SphericalTransformation(pdom, hdom, self._space) self.lmax = hspc.lmax
self.mmax = hspc.mmax
self.sjob = sharpjob_d()
self.sjob.set_triangular_alm_info(self.lmax, self.mmax)
if isinstance(pspc, GLSpace):
self.sjob.set_Gauss_geometry(pspc.nlat, pspc.nlon)
else:
self.sjob.set_Healpix_geometry(pspc.nside)
if adom.harmonic:
self._capability = self.TIMES | self.ADJOINT_TIMES
else:
self._capability = (self.INVERSE_TIMES |
self.INVERSE_ADJOINT_TIMES)
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
if np.issubdtype(x.dtype, np.complexfloating): if np.issubdtype(x.dtype, np.complexfloating):
res = (self._trafo.transform(x.real) + return (self._applyfunc(x.real, mode) +
1j * self._trafo.transform(x.imag)) 1j*self._applyfunc(x.imag, mode))
else:
return self._applyfunc(x, mode)
def _apply_cartesian(self, x, mode):
"""
RG -> RG transform method.
Parameters
----------
x : Field
The field to be transformed
"""
from pyfftw.interfaces.numpy_fft import fftn
axes = x.domain.axes[self._space]
tdom = self._target if x.domain == self._domain else self._domain
oldax = dobj.distaxis(x.val)
if oldax not in axes: # straightforward, no redistribution needed
ldat = dobj.local_data(x.val)
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 = 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: else:
res = self._trafo.transform(x) fct = self._target[self._space].scalar_dvol()
return res if fct != 1:
Tval *= fct
return Tval
def _slice_p2h(self, inp):
rr = self.sjob.alm2map_adjoint(inp)
assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax)
res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype)
res[0:self.lmax+1] = rr[0:self.lmax+1].real
res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real
res[self.lmax+2::2] = np.sqrt(2)*rr[self.lmax+1:].imag
return res/np.sqrt(np.pi*4)
def _slice_h2p(self, inp):
res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype)
assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax)
res[0:self.lmax+1] = inp[0:self.lmax+1]
res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] +
1j*inp[self.lmax+2::2])
res = self.sjob.alm2map(res)
return res/np.sqrt(np.pi*4)
def _apply_spherical(self, x, mode):
axes = x.domain.axes[self._space]
axis = axes[0]
tval = x.val
if dobj.distaxis(tval) == axis:
tval = dobj.redistribute(tval, nodist=(axis,))
distaxis = dobj.distaxis(tval)
p2h = not x.domain[self._space].harmonic
tdom = self._target if x.domain == self._domain else self._domain
func = self._slice_p2h if p2h else self._slice_h2p
idat = dobj.local_data(tval)
odat = np.empty(dobj.local_shape(tdom.shape, distaxis=distaxis),
dtype=x.dtype)
for slice in utilities.get_slice_list(idat.shape, axes):
odat[slice] = func(idat[slice])
odat = dobj.from_local_data(tdom.shape, odat, distaxis)
if distaxis != dobj.distaxis(x.val):
odat = dobj.redistribute(odat, dist=dobj.distaxis(x.val))
return Field(tdom, odat)
@property @property
def domain(self): def domain(self):
...@@ -103,7 +222,4 @@ class FFTOperator(LinearOperator): ...@@ -103,7 +222,4 @@ class FFTOperator(LinearOperator):
@property @property
def capability(self): def capability(self):
res = self.TIMES | self.ADJOINT_TIMES return self._capability
if self._trafo.unitary:
res |= self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES
return res
# 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-2017 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 division
import numpy as np
from .. import utilities
from .. import dobj
from ..field import Field
from ..spaces.gl_space import GLSpace
class Transformation(object):
def __init__(self, pdom, hdom, space):
self.pdom = pdom
self.hdom = hdom
self.space = space
class RGRGTransformation(Transformation):
def __init__(self, pdom, hdom, space):
import pyfftw
super(RGRGTransformation, self).__init__(pdom, hdom, space)
pyfftw.interfaces.cache.enable()
# correct for forward/inverse fft.
# naively one would set power to 0.5 here in order to
# apply effectively a factor of 1/sqrt(N) to the field.
# BUT: the pixel volumes of the domain and codomain are different.
# Hence, in order to produce the same scalar product, power==1.
self.fct_p2h = pdom[space].scalar_dvol()
self.fct_h2p = 1./(pdom[space].scalar_dvol()*hdom[space].dim)
@property
def unitary(self):
return True
def transform(self, x):
"""
RG -> RG transform method.
Parameters
----------
x : Field
The field to be transformed
"""
from pyfftw.interfaces.numpy_fft import fftn
axes = x.domain.axes[self.space]
p2h = x.domain == self.pdom
tdom = self.hdom if p2h else self.pdom
oldax = dobj.distaxis(x.val)
if oldax not in axes: # straightforward, no redistribution needed
ldat = dobj.local_data(x.val)
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
if True:
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)
# new, experimental code
if True:
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 = 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)
else:
tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=oldax)
tmp = dobj.redistribute(tmp, nodist=(oldax,))
newax = dobj.distaxis(tmp)
ldat = dobj.local_data(tmp)
ldat = fftn(ldat, axes=(oldax,))
ldat = ldat.real+ldat.imag
tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax)
tmp = dobj.redistribute(tmp, dist=oldax)
else:
tmp = dobj.redistribute(x.val, nodist=(oldax,))
newax = dobj.distaxis(tmp)
ldat = dobj.local_data(tmp)
ldat = fftn(ldat, axes=(oldax,))
tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax)
tmp = dobj.redistribute(tmp, dist=oldax)
rem_axes = tuple(i for i in axes if i != oldax)
ldat = dobj.local_data(tmp)
ldat = fftn(ldat, axes=rem_axes)
ldat = ldat.real+ldat.imag
tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=oldax)
Tval = Field(tdom, tmp)
fct = self.fct_p2h if p2h else self.fct_h2p
if fct != 1:
Tval *= fct
return Tval
class SphericalTransformation(Transformation):
def __init__(self, pdom, hdom, space):
super(SphericalTransformation, self).__init__(pdom, hdom, space)
from pyHealpix import sharpjob_d