Commit 5658335b authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'NIFTy_2' into 'master'

Nifty 2

Merging the changes from theos into master. 

See merge request !3
parents 35b4c239 b41e5ad8
*.html
*.c
*.o
*.so
*.pyc
......@@ -10,26 +11,5 @@
.git/
.svn/
.spyderproject
.document
build
rg/*
!rg/__init__.py
!rg/nifty_rg.py
!rg/nifty_power_conversion_rg.py
!rg/gfft_rg.py
!rg/powerspectrum.py
lm/*
!lm/__init__.py
!lm/nifty_lm.py
!lm/nifty_power_conversion_lm.py
demos/*
!demos/__init__.py
!demos/demos_core.py
!demos/demo_faraday.py
!demos/demo_faraday_map.npy
!demos/demo_excaliwir.py
!demos/demo_wf1.py
!demos/demo_wf2.py
!demos/demo_wf3.py
\ No newline at end of file
# 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
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
......@@ -20,24 +20,66 @@
## along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
from nifty_core import * ## imports `about`
from nifty_cmaps import *
from nifty_power import *
from nifty_tools import *
from nifty_explicit import *
import matplotlib as mpl
mpl.use('Agg')
import dummys
from keepers import about,\
global_dependency_injector,\
global_configuration
from nifty_cmaps import ncmap
from nifty_core import space,\
point_space,\
field
from nifty_mpi_data import distributed_data_object, d2o_librarian
from nifty_random import random
from nifty_simple_math import *
from nifty_utilities import *
from nifty_paradict import space_paradict,\
point_space_paradict,\
nested_space_paradict
from operators import *
## optional submodule `rg`
try:
from rg import *
from rg import rg_space,\
power_backward_conversion_rg,\
power_forward_conversion_rg
from nifty_paradict import rg_space_paradict
except(ImportError):
pass
## optional submodule `lm`
try:
from lm import *
from lm import lm_space,\
power_backward_conversion_lm,\
power_forward_conversion_lm
from nifty_paradict import lm_space_paradict
try:
from lm import gl_space
from nifty_paradict import gl_space_paradict
except(ImportError):
pass
try:
from lm import hp_space
from nifty_paradict import hp_space_paradict
except(ImportError):
pass
except(ImportError):
pass
from demos import *
from pickling import *
from demos import get_demo_dir
from pickling import _pickle_method, _unpickle_method
#import pyximport; pyximport.install(pyimport = True)
......@@ -19,5 +19,5 @@
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
from demos_core import *
from demos_core import get_demo_dir
......@@ -32,6 +32,13 @@
"""
from __future__ import division
import matplotlib as mpl
mpl.use('Agg')
import imp
#nifty = imp.load_module('nifty', None,
# '/home/steininger/Downloads/nifty', ('','',5))
from nifty import *
......@@ -39,7 +46,7 @@ from nifty import *
class problem(object):
def __init__(self, x_space, s2n=2, **kwargs):
def __init__(self, x_space, s2n=6, **kwargs):
"""
Sets up a Wiener filter problem.
......@@ -51,14 +58,16 @@ class problem(object):
Signal-to-noise ratio (default: 2).
"""
self.store = []
## set signal space
self.z = x_space
## set conjugate space
self.k = self.z.get_codomain()
self.k.set_power_indices(**kwargs)
#self.k.power_indices.set_default()
#self.k.set_power_indices(**kwargs)
## set some power spectrum
self.power = (lambda k: 42 / (k + 1) ** 3)
self.power = (lambda k: 42 / (k + 1) ** 2)
## define signal covariance
self.S = power_operator(self.k, spec=self.power, bare=True)
......@@ -73,7 +82,8 @@ class problem(object):
d_space = self.R.target
## define noise covariance
self.N = diagonal_operator(d_space, diag=abs(s2n) * self.s.var(), bare=True)
#self.N = diagonal_operator(d_space, diag=abs(s2n) * self.s.var(), bare=True)
self.N = diagonal_operator(d_space, diag=abs(s2n), bare=True)
## define (plain) projector
self.Nj = projection_operator(d_space)
## generate noise
......@@ -83,9 +93,12 @@ class problem(object):
self.d = self.R(self.s) + n
## define information source
self.j = self.R.adjoint_times(self.N.inverse_times(self.d), target=self.k)
#self.j = self.R.adjoint_times(self.N.inverse_times(self.d), target=self.k)
self.j = self.R.adjoint_times(self.N.inverse_times(self.d))
## define information propagator
self.D = propagator_operator(S=self.S, N=self.N, R=self.R)
self.D = propagator_operator(S=self.S,
N=self.N,
R=self.R)
## reserve map
self.m = None
......@@ -150,32 +163,74 @@ class problem(object):
## pre-compute denominator
denominator = self.k.power_indices["rho"] + 2 * (alpha - 1 + abs(epsilon))
self.save_signal_and_data()
## iterate
i = 0
iterating = True
while(iterating):
## reconstruct map
self.m = self.D(self.j, W=self.S, tol=1E-3, note=False)
self.m = self.D(self.j, W=self.S, tol=1E-3, note=True)
if(self.m is None):
break
#print'Reconstructed m'
## reconstruct power spectrum
tr_B1 = self.Sk.pseudo_tr(self.m) ## == Sk(m).pseudo_dot(m)
print 'Calculated trace B1'
print ('tr_b1', tr_B1)
tr_B2 = self.Sk.pseudo_tr(self.D, loop=True)
numerator = 2 * q + tr_B1 + abs(delta) * tr_B2 ## non-bare(!)
print 'Calculated trace B2'
print ('tr_B2', tr_B2)
numerator = 2 * q + tr_B1 + tr_B2 * abs(delta) ## non-bare(!)
power = numerator / denominator
print ('numerator', numerator)
print ('denominator', denominator)
print ('power', power)
print 'Calculated power'
#power = np.clip(power, 0.00000001, np.max(power))
self.store += [{'tr_B1': tr_B1,
'tr_B2': tr_B2,
'num': numerator,
'denom': denominator}]
## check convergence
dtau = log(power / self.S.get_power(), base=self.S.get_power())
print ('dtau', np.max(np.abs(dtau)))
iterating = (np.max(np.abs(dtau)) > 2E-2)
print max(np.abs(dtau))
#printmax(np.abs(dtau))
self.save_map(i)
i += 1
## update signal covariance
self.S.set_power(power, bare=False) ## auto-updates D
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def save_signal_and_data(self):
self.s.plot(title="signal", save="img/signal.png")
try:
d_ = field(self.z, val=self.d.val, target=self.k)
d_.plot(title="data", vmin=self.s.min(), vmax=self.s.max(),
save="img/data.png")
except:
pass
def save_map(self, index=0):
# save map
if(self.m is None):
pass
else:
self.m.plot(title="reconstructed map",
vmin=self.s.min(), vmax=self.s.max(),
save="img/map_"+str(index)+".png")
self.m.plot(power=True, mono=False,
other=(self.power, self.S.get_power()),
nbin=None, binbounds=None, log=False,
save='img/map_power_'+str(index)+".png")
def plot(self):
"""
Produces plots.
......@@ -199,36 +254,39 @@ class problem(object):
##=============================================================================
##-----------------------------------------------------------------------------
#
if(__name__=="__main__"):
# pl.close("all")
## define signal space
x_space = rg_space(128)
## setup problem
p = problem(x_space, log=True)
## solve problem given some power spectrum
p.solve()
## solve problem
p.solve_critical()
p.plot()
## retrieve objects
k_space = p.k
power = p.power
S = p.S
Sk = p.Sk
s = p.s
R = p.R
d_space = p.R.target
N = p.N
Nj = p.Nj
d = p.d
j = p.j
D = p.D
m = p.m
x = rg_space((1280), zerocenter=True)
p = problem(x, log = False)
about.warnings.off()
## pl.close("all")
#
# ## define signal space
# x_space = rg_space(128)
#
# ## setup problem
# p = problem(x_space, log=True)
# ## solve problem given some power spectrum
# p.solve()
# ## solve problem
# p.solve_critical()
#
# p.plot()
#
# ## retrieve objects
# k_space = p.k
# power = p.power
# S = p.S
# Sk = p.Sk
# s = p.s
# R = p.R
# d_space = p.R.target
# N = p.N
# Nj = p.Nj
# d = p.d
# j = p.j
# D = p.D
# m = p.m
##-----------------------------------------------------------------------------
......@@ -101,7 +101,7 @@ def run(projection=False, power=False):
z_space = rg_space([768, 384], dist=[1/360, 1/180])
m5 = m1.transform(y_space)
m5.cast_domain(z_space)
m5.set_val(np.roll(m5.val[::-1, ::-1], y_space.nlon()//2, axis=1)) # rearrange value array
m5.set_val(np.roll(m5.val[::-1, ::-1], y_space.paradict['nlon']//2, axis=1)) # rearrange value array
m5.plot(title=r"$m$ on a regular 2D grid", **nicely)
if(power):
......
# -*- coding: utf-8 -*-
from nifty import *
if __name__ == "__main__":
about.warnings.off()
shape = (256, 256)
x_space = rg_space(shape)
k_space = x_space.get_codomain()
power = lambda k: 42/((1+k*shape[0])**3)
S = power_operator(k_space, codomain=x_space, spec=power)
s = S.get_random_field(domain=x_space)
#n_points = 360.
#starts = [[(np.cos(i/n_points*np.pi)+1)*shape[0]/2.,
# (np.sin(i/n_points*np.pi)+1)*shape[0]/2.] for i in xrange(int(n_points))]
#starts = list(np.array(starts).T)
#
#ends = [[(np.cos(i/n_points*np.pi + np.pi)+1)*shape[0]/2.,
# (np.sin(i/n_points*np.pi + np.pi)+1)*shape[0]/2.] for i in xrange(int(n_points))]
#ends = list(np.array(ends).T)
def make_los(n=10, angle=0, d=1):
starts_list = []
ends_list = []
for i in xrange(n):
starts_list += [[(-0.2)*d, (-0.2 + 1.2*i/n)*d]]
ends_list += [[(1.2)*d, (-0.2 + 1.2*i/n)*d]]
starts_list = np.array(starts_list)
ends_list = np.array(ends_list)
rot_matrix = np.array([[np.cos(angle), -np.sin(angle)],
[np.sin(angle), np.cos(angle)]])
starts_list = rot_matrix.dot(starts_list.T-0.5*d).T+0.5*d
ends_list = rot_matrix.dot(ends_list.T-0.5*d).T+0.5*d
return (starts_list, ends_list)
temp_coords = (np.empty((0, 2)), np.empty((0, 2)))
n = 250
m = 250
for alpha in [np.pi/n*j for j in xrange(n)]:
temp = make_los(n=m, angle=alpha)
temp_coords = np.concatenate([temp_coords, temp], axis=1)
starts = list(temp_coords[0].T)
ends = list(temp_coords[1].T)
R = los_response(x_space, starts=starts, ends=ends,
sigmas_up=0.1, sigmas_low=0.1)
d_space = R.target
N = diagonal_operator(d_space, diag=s.var()/100000, bare=True)
n = N.get_random_field(domain=d_space)
d = R(s) + n
j = R.adjoint_times(N.inverse_times(d))
D = propagator_operator(S=S, N=N, R=R)
m = D(j, W=S, tol=1E-14, limii=100, note=True)
s.plot(title="signal", save='1_plot_s.png')
s.plot(save='plot_s_power.png', power=True, other=power)
j.plot(save='plot_j.png')
#d_ = field(x_space, val=d.val, target=k_space)
#d_.plot(title="data", vmin=s.min(), vmax=s.max(), save='plot_d.png')
m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save='1_plot_m.png')
m.plot(title="reconstructed map", save='2_plot_m.png')
m.plot(save='plot_m_power.png', power=True, other=power)
\ No newline at end of file
......@@ -32,39 +32,86 @@
"""
from __future__ import division
import matplotlib as mpl
mpl.use('Agg')
import gc
#import imp
#nifty = imp.load_module('nifty', None,
# '/home/steininger/Downloads/nifty', ('','',5))
from nifty import * # version 0.8.0
if __name__ == "__main__":