Commit 2f17dc3e authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into 'spherical_plots'

Master

See merge request !127
parents 7d4f9288 046142bf
Pipeline #12555 passed with stage
in 6 minutes and 23 seconds
......@@ -78,11 +78,6 @@ The current version of Nifty3 can be obtained by cloning the repository:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
and switching to the "master" branch:
cd NIFTy
git checkout master
### Installation on Ubuntu
This is for you if you want to install NIFTy on your personal computer
......@@ -98,14 +93,13 @@ Starting with a fresh Ubuntu installation move to a folder like
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure --prefix=$HOME/.local && make -j4 && make install
autoreconf -i && ./configure --prefix=$HOME/.local --enable-openmp --enable-native-optimizations && make -j4 && make install
cd ..
- Finally, NIFTy:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
cd NIFTy
git checkout master
python setup.py install --user
cd ..
......@@ -135,7 +129,7 @@ may cause trouble.
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure --prefix=`python-config --prefix` && make -j4 && sudo make install
autoreconf -i && ./configure --prefix=`python-config --prefix` --enable-openmp --enable-native-optimizations && make -j4 && sudo make install
cd ..
(The third command installs the package system-wide. User-specific
......@@ -146,7 +140,6 @@ may cause trouble.
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
cd NIFTy
git checkout master
python setup.py install --user
cd ..
......@@ -197,4 +190,4 @@ The NIFTY package is licensed under the
[1] Selig et al., "NIFTY - Numerical Information Field Theory - a
versatile Python library for signal inference", [A&A, vol. 554, id.
A26](http://dx.doi.org/10.1051/0004-6361/201321236), 2013;
[arXiv:1301.4499](http://www.arxiv.org/abs/1301.4499)
\ No newline at end of file
[arXiv:1301.4499](http://www.arxiv.org/abs/1301.4499)
#!/bin/bash
wget https://api.github.com/repos/h5py/h5py/tags -O - | grep tarball_url | grep -v rc | head -n 1 | cut -d '"' -f 4 | wget -i - -O h5py.tar.gz
tar xzf h5py.tar.gz
cd h5py-h5py*
export CC=mpicc
export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi
python setup.py configure --mpi
python setup.py build
python setup.py install
cd ..
rm -r h5py-h5py*
rm h5py.tar.gz
export HDF5_MPI="ON"
pip install --no-binary=h5py h5py
numpy
cython
mpi4py
matplotlib
plotly
......
numpy
cython
nose
parameterized
coverage
......
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().real)], 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(32)
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
Fft_op = FFTOperator(domain=x1, target=k1,
domain_dtype=np.float64,
target_dtype=np.complex128)
R_op = ResponseOperator(x1, sigma=[length_convolution],
exposure=[exposure])
# drawing a random field
sk = p_field.power_synthesize(real_signal=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
j = Fft_op.times(R_op.adjoint_times(N_op.inverse_times(d)))
D = HarmonicPropagatorOperator(S=S_op, N=N_op, R=R_op)
mk = D(j)
m = Fft_op.adjoint_times(mk)
# z={}
# z["signal"] = s
# z["reconstructed_map"] = m
# z["data"] = d
# z["lambda"] = R_op(s)
# z["j"] = j
#
# plot_maps(z, "Wiener_filter.html")
......@@ -17,6 +17,8 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
import itertools
import numpy as np
from keepers import Versionable,\
......@@ -104,6 +106,7 @@ class Field(Loggable, Versionable, object):
distributed_data_object
"""
# ---Initialization methods---
def __init__(self, domain=None, val=None, dtype=None,
......@@ -358,7 +361,7 @@ class Field(Loggable, Versionable, object):
self.domain_axes[space_index])
harmonic_domain = self.domain[space_index]
power_domain = PowerSpace(harmonic_domain=harmonic_domain,
power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
......@@ -512,7 +515,7 @@ class Field(Loggable, Versionable, object):
result_domain = list(self.domain)
for power_space_index in spaces:
power_space = self.domain[power_space_index]
harmonic_domain = power_space.harmonic_domain
harmonic_domain = power_space.harmonic_partner
result_domain[power_space_index] = harmonic_domain
# create random samples: one or two, depending on whether the
......@@ -554,14 +557,12 @@ class Field(Loggable, Versionable, object):
inplace=True)
if real_signal:
for power_space_index in spaces:
harmonic_domain = result_domain[power_space_index]
result_val_list = [harmonic_domain.hermitian_decomposition(
result_val,
axes=result.domain_axes[power_space_index],
preserve_gaussian_variance=True)[0]
for (result, result_val)
in zip(result_list, result_val_list)]
result_val_list = [self._hermitian_decomposition(
result_domain,
result_val,
spaces,
result_list[0].domain_axes)[0]
for result_val in result_val_list]
# store the result into the fields
[x.set_val(new_val=y, copy=False) for x, y in
......@@ -574,6 +575,46 @@ class Field(Loggable, Versionable, object):
return result
@staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition(
val,
domain_axes[spaces[0]])
# hermitianize all remaining spaces using the iterative formula
for space in xrange(1, len(spaces)):
(hh, ha) = \
domain[space].hermitian_decomposition(h, domain_axes[space])
(ah, aa) = \
domain[space].hermitian_decomposition(a, domain_axes[space])
c = (hh - ha - ah + aa).conjugate()
h = (val + c)/2.
a = (val - c)/2.
# correct variance
fixed_points = [domain[i].hermitian_fixed_points() for i in spaces]
# check if there was at least one flipping during hermitianization
flipped_Q = np.any([fp is not None for fp in fixed_points])
# if the array got flipped, correct the variance
if flipped_Q:
h *= np.sqrt(2)
a *= np.sqrt(2)
fixed_points = [[fp] if fp is None else fp for fp in fixed_points]
for product_point in itertools.product(*fixed_points):
slice_object = np.array((slice(None), )*len(val.shape),
dtype=np.object)
for i, sp in enumerate(spaces):
point_component = product_point[i]
if point_component is None:
point_component = slice(None)
slice_object[list(domain_axes[sp])] = point_component
slice_object = tuple(slice_object)
h[slice_object] /= np.sqrt(2)
a[slice_object] /= np.sqrt(2)
return (h, a)
def _spec_to_rescaler(self, spec, result_list, power_space_index):
power_space = self.domain[power_space_index]
......
......@@ -94,9 +94,12 @@ class LMSpace(Space):
hermitian_part = x.copy_empty()
anti_hermitian_part = x.copy_empty()
hermitian_part[:] = x.real
anti_hermitian_part[:] = x.imag
anti_hermitian_part[:] = x.imag * 1j
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self):
return None
# ---Mandatory properties and methods---
@property
......
......@@ -114,6 +114,11 @@ class PowerSpace(Space):
self._pundex = power_index['pundex']
self._k_array = power_index['k_array']
if self.config['nbin'] is not None:
if self.config['nbin'] > len(self.kindex):
self.logger.warn("nbin was set to a value being larger than "
"the length of kindex!")
def pre_cast(self, x, axes):
""" Casts power spectrum functions to discretized power spectra.
......
......@@ -110,7 +110,7 @@ class RGSpace(Space):
hermitian_part /= 2.
# use subtraction since it is faster than flipping another time
anti_hermitian_part = (x-hermitian_part)/1j
anti_hermitian_part = (x-hermitian_part)
if preserve_gaussian_variance:
hermitian_part, anti_hermitian_part = \
......@@ -120,6 +120,18 @@ class RGSpace(Space):
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self):
shape = self.shape
mid_index = np.array(shape)//2
ndlist = [2 if (shape[i] % 2 == 0) else 1 for i in xrange(len(shape))]
ndlist = tuple(ndlist)
odd_axes_list = np.array([1 if (shape[i] % 2 == 1) else 0
for i in xrange(len(shape))])
fixed_points = []
for i in np.ndindex(ndlist):
fixed_points += [tuple((i+odd_axes_list) * mid_index)]
return fixed_points
def _hermitianize_correct_variance(self, hermitian_part,
anti_hermitian_part, axes):
# Correct the variance by multiplying sqrt(2)
......@@ -145,8 +157,9 @@ class RGSpace(Space):
return hermitian_part, anti_hermitian_part
def _hermitianize_inverter(self, x, axes):
shape = x.shape
# calculate the number of dimensions the input array has
dimensions = len(x.shape)
dimensions = len(shape)
# prepare the slicing object which will be used for mirroring
slice_primitive = [slice(None), ] * dimensions
# copy the input data
......@@ -158,11 +171,17 @@ class RGSpace(Space):
# flip in the desired directions
for i in axes:
slice_picker = slice_primitive[:]
slice_picker[i] = slice(1, None, None)
if shape[i] % 2 == 0:
slice_picker[i] = slice(1, None, None)
else:
slice_picker[i] = slice(None)
slice_picker = tuple(slice_picker)
slice_inverter = slice_primitive[:]
slice_inverter[i] = slice(None, 0, -1)
if shape[i] % 2 == 0:
slice_inverter[i] = slice(None, 0, -1)
else:
slice_inverter[i] = slice(None, None, -1)
slice_inverter = tuple(slice_inverter)
try:
......
......@@ -21,31 +21,18 @@ import unittest
from numpy.testing import assert_equal
from keepers import Repository
from test.common import expand, generate_spaces
from nose.plugins.skip import SkipTest
import os
try:
import h5py
except ImportError:
h5py_available = False
else:
h5py_available = True
if h5py_available:
class SpaceSerializationTests(unittest.TestCase):
# variable to hold the repository
_repo = None
@classmethod
def setUpClass(cls):
cls._repo = Repository('test.h5')
@expand([[space] for space in generate_spaces()])
def test_serialization(self, space):
self._repo.add(space, 'space')
self._repo.commit()
assert_equal(space, self._repo.get('space'))
@classmethod
def tearDownClass(cls):
os.remove('test.h5')
class SpaceSerializationTests(unittest.TestCase):
@expand([[space] for space in generate_spaces()])
def test_serialization(self, space):
try:
import h5py
except ImportError:
raise SkipTest
repo = Repository('test.h5')
repo.add(space, 'space')
repo.commit()
assert_equal(space, repo.get('space'))
os.remove('test.h5')
......@@ -116,7 +116,7 @@ class LMSpaceFunctionalityTests(unittest.TestCase):
real)
assert_almost_equal(
l.hermitian_decomposition(distributed_data_object(x))[1],
imag)
imag*1j)
@expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected):
......
......@@ -162,7 +162,7 @@ def get_hermitian_configs():
[0.37928780+0.j, 0.33013206+0.26511434j, 0.48235714+0.j,
0.33013206-0.26511434j],
[0.47773437+0.17342852j, 0.31059374+0.02379381j, 0.64003551-0.23838932j,
0.27824437-0.0197826j]])
0.27824437-0.0197826j]])*1j
h_1_x = np.array([
[[0.23987021+0.41617749j, 0.34605012+0.55462234j, 0.07947035+0.73360723j,
......@@ -205,7 +205,7 @@ def get_hermitian_configs():
0.64466122-0.10318736j],
[0.19076379+0.j, 0.49389371+0.03664104j, 0.85091112+0.j,
0.49389371-0.03664104j]]
])
])*1j
return [
[h_0_x, None, h_0_res_real, h_0_res_imag],
[h_1_x, (2,), h_1_res_real, h_1_res_imag]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment