Merge branch 'master' into field_docu

.gitlab-ci.yml

@@ -13,7 +13,7 @@ before_script:
   - apt-get update
   - >
     apt-get install -y build-essential python python-pip python-dev git
-    autoconf gsl-bin libgsl-dev wget python-numpy cython
+    autoconf gsl-bin libgsl-dev wget python-numpy
   - pip install --upgrade -r ci/requirements_base.txt
   - chmod +x ci/*.sh
 

Makefile

-# This Makefile implements common tasks needed by developers
-# A list of implemented rules can be obtained by the command "make help"
-
-.DEFAULT_GOAL=build
-.PHONY .SILENT : help
-help :
-	echo
-	echo "    Implemented targets:"
-	echo
-	echo "    build        build pypmc for python2 and python3"
-	echo "    buildX       build pypmc for pythonX only where X is one of {2,3}"
-	echo "    build-sdist  build pypmc from the dist directory (python 2 and 3)"
-	echo "    build-sdistX build pypmc from the dist directory (pythonX, X in {2,3})"
-	echo "    check        use nosetests to test pypmc with python 2.7 and 3"
-	echo "    checkX       use nosetests to test pypmc with python 2.7 or 3,"
-	echo "                 where X is one of {2,3}"
-	echo "    check-fast   use nosetests to run only quick tests of pypmc"
-	echo "                 using nosetests-2.7 and nosetests3"
-	echo "    check-sdist  use nosetests-2.7 and nosetests3 to test the distribution"
-	echo "                 generated by 'make sdist'"
-	echo "    check-sdistX use nosetests-2.7 or nosetests3 to test the distribution"
-	echo "                 generated by 'make sdist', where X is one of {2,3}"
-	echo "    clean        delete compiled and temporary files"
-	echo "    coverage     produce and show a code coverage report"
-	echo "                 Note: Cython modules cannot be analyzed"
-	echo "    distcheck    runs 'check', check-sdist', 'run-examples' and"
-	echo "                 opens a browser with the built documentation"
-	echo "    doc          build the html documentation using sphinx"
-	echo "    doc-pdf      build the pdf documentation using sphinx"
-	echo "    help         show this message"
-	echo "    run-examples run all examples using python 2 and 3"
-	echo "    sdist        make a source distribution"
-	echo "    show-todos   show todo marks in the source code"
-	echo
-
-.PHONY : clean
-clean:
-	#remove build doc
-	rm -rf ./doc/_build
-
-	#remove .pyc files created by python 2.7
-	rm -f ./*.pyc
-	find -P . -name '*.pyc' -delete
-
-	#remove .pyc files crated by python 3
-	rm -rf ./__pycache__
-	find -P . -name __pycache__ -delete
-
-	#remove build folder in root directory
-	rm -rf ./build
-
-	#remove cythonized C source and object files
-	find -P . -name '*.c' -delete
-
-	#remove variational binaries only if command line argument specified
-	find -P . -name '*.so' -delete
-
-	#remove backup files
-	find -P . -name '*~' -delete
-
-	#remove files created by coverage
-	rm -f .coverage
-	rm -rf coverage
-
-	# remove egg info
-	rm -rf pypmc.egg-info
-
-	# remove downloaded seutptools
-	rm -f setuptools-3.3.zip
-
-	# remove dist/
-	rm -rf dist
-
-.PHONY : build
-build : build2
-
-.PHONY : build2
-build2 :
-	python2 setup.py build_ext --inplace
-
-.PHONY :
-check : check2
-
-.PHONY : check2
-check2 : build2
-	@ # run tests
-	nosetests-2.7 --processes=-1 --process-timeout=60
-
-	# run tests in parallel
-	mpirun -n 2 nosetests-2.7
-
-.PHONY : check-fast
-check-fast : build
-	nosetests-2.7 -a '!slow' --processes=-1 --process-timeout=60
-	nosetests3    -a '!slow' --processes=-1 --process-timeout=60
-
-.PHONY : .build-system-default
-.build-system-default :
-	python setup.py build_ext --inplace
-
-.PHONY : doc
-doc : .build-system-default
-	cd doc && make html
-
-.PHONY : doc-pdf
-doc-pdf : .build-system-default
-	cd doc; make latexpdf
-
-.PHONY : run-examples
-run-examples : build
-	cd examples ; \
-	for file in $$(ls) ; do \
-	    echo running $${file} with python2 && \
-	    python2 $${file} && \
-	    echo running $${file} with python3 && \
-	    python3 $${file} && \
-	    \
-	    # execute with mpirun if mpi4py appears in the file \
-	    if grep -Fq 'mpi4py' $${file} ; then \
-		echo "$${file}" is mpi parallelized && \
-		echo running $${file} in parallel with python2 && \
-		mpirun -n 2 python2 $${file} && \
-		echo running $${file} in parallel with python3 && \
-		mpirun -n 2 python3 $${file} ; \
-	    fi \
-	; \
-	done
-
-.PHONY : sdist
-sdist :
-	python setup.py sdist
-
-.PHONY : build-sdist
-build-sdist : build-sdist2 build-sdist3
-
-./dist/pypmc*/NUL : sdist
-	cd dist && tar xaf *.tar.gz && cd *
-
-.PHONY : build-sdist2
-build-sdist2 : ./dist/pypmc*/NUL
-	cd dist/pypmc* && python2 setup.py build
-
-.PHONY : build-sdist3
-build-sdist3 : ./dist/pypmc*/NUL
-	cd dist/pypmc* && python3 setup.py build
-
-.PHONY : check-sdist
-check-sdist : check-sdist2 check-sdist3
-
-.PHONY : check-sdist2
-check-sdist2 : build-sdist2
-	cd dist/*/build/lib*2.7 && \
-	nosetests-2.7 --processes=-1 --process-timeout=60 && \
-	mpirun -n 2 nosetests-2.7
-
-.PHONY : check-sdist3
-check-sdist3 : build-sdist3
-	cd dist/*/build/lib*3.* && \
-	nosetests3 --processes=-1 --process-timeout=60 && \
-	mpirun -n 2 nosetests3
-
-.PHONY : distcheck
-distcheck : check check-sdist doc
-	@ # execute "run-examples" after all other recipes makes are done
-	make run-examples
-	xdg-open link_to_documentation
-
-.PHONY : show-todos
-grep_cmd = ack-grep -i --no-html --no-cc [^"au""sphinx.ext."]todo
-show-todos :
-	@ # suppress errors here
-	@ # note that no todo found is considered as error
-	$(grep_cmd) doc ; \
-	$(grep_cmd) pypmc ; \
-	$(grep_cmd) examples ; echo \
-
-.PHONY : coverage
-coverage : .build-system-default
-	rm -rf coverage
-	nosetests --with-coverage --cover-package=nifty --cover-html --cover-html-dir=coverage
-	xdg-open coverage/index.html

PKG-INFO

-Metadata-Version: 1.0
-Name: ift_nifty
-Version: 1.0.6
-Summary: Numerical Information Field Theory
-Home-page: http://www.mpa-garching.mpg.de/ift/nifty/
-Author: Theo Steininger
-Author-email: theos@mpa-garching.mpg.de
-License: GPLv3
-Description: UNKNOWN
-Platform: UNKNOWN

README.md

@@ -15,7 +15,7 @@ Summary
 a versatile library designed to enable the development of signal
 inference algorithms that operate regardless of the underlying spatial
 grid and its resolution. Its object-oriented framework is written in
-Python, although it accesses libraries written in Cython, C++, and C for
+Python, although it accesses libraries written in C++ and C for
 efficiency.
 
 NIFTY offers a toolkit that abstracts discretized representations of
@@ -71,7 +71,6 @@ Installation
 
 -   [Python](http://www.python.org/) (v2.7.x)
     -   [NumPy](http://www.numpy.org/)
-    -   [Cython](http://cython.org/)
 
 ### Download
 
@@ -95,7 +94,7 @@ Starting with a fresh Ubuntu installation move to a folder like
 
 -   Using pip install numpy etc...:
 
-        sudo pip install numpy cython
+        sudo pip install numpy
 
 -   Install pyHealpix:
 
@@ -147,10 +146,9 @@ MacPorts, missing ones need to be installed manually. It may also be
 mentioned that one should only use one package manager, as multiple ones
 may cause trouble.
 
--   Install basic packages numpy and cython:
+-   Install numpy:
 
         sudo port install py27-numpy
-        sudo port install py27-cython
 
 -   Install pyHealpix:
 

ci/requirements.txt

 numpy
-cython
 mpi4py
 matplotlib
 plotly

demos/wiener_filter.py

@@ -11,16 +11,36 @@ rank = comm.rank
 if __name__ == "__main__":
 
     distribution_strategy = 'not'
-
+    
+    #Setting up physical constants
+    #total length of Interval or Volume the field lives on, e.g. in meters
+    L = 2.
+    #typical distance over which the field is correlated (in same unit as L)
+    correlation_length = 0.1
+    #variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
+    field_variance = 2.
+    #smoothing length that response (in same unit as L)
+    response_sigma = 0.1
+    
+    #defining resolution (pixels per dimension)
+    N_pixels = 512
+    
+    #Setting up derived constants
+    k_0 = 1./correlation_length
+    #note that field_variance**2 = a*k_0/4. for this analytic form of power
+    #spectrum
+    a = field_variance**2/k_0*4.
+    pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
+    pixel_width = L/N_pixels
+    
     # Setting up the geometry
-    s_space = RGSpace([512, 512])
+    s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
     fft = FFTOperator(s_space)
     h_space = fft.target[0]
     p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
 
 
     # Creating the mock data
-    pow_spec = (lambda k: 42 / (k + 1) ** 3)
 
     S = create_power_operator(h_space, power_spectrum=pow_spec,
                               distribution_strategy=distribution_strategy)
@@ -30,7 +50,7 @@ if __name__ == "__main__":
     sh = sp.power_synthesize(real_signal=True)
     ss = fft.inverse_times(sh)
 
-    R = SmoothingOperator(s_space, sigma=0.1)
+    R = SmoothingOperator(s_space, sigma=response_sigma)
 
     signal_to_noise = 1
     N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)

demos/wiener_filter_harmonic - Kopie.py

+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")
+
+

demos/wiener_filter_harmonic.py

+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())], 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(64)
+    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
+    # adjust dtype_target probperly
+    Fft_op = FFTOperator(domain=x1, target=k1,
+                        domain_dtype=np.float64,
+                        target_dtype=np.float64)
+    R_op = ResponseOperator(x1, sigma=[length_convolution],
+                            exposure=[exposure])
+
+    # drawing a random field
+    sk = p_field.power_synthesize(real_power=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)
+    #
+    # plot_maps(z, "Wiener_filter.html")
+
+

demos/wiener_filter_unit.py

+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())], 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