Skip to content

GitLab

Commit 69a98326 authored by Theo Steininger's avatar Theo Steininger
Browse files
Options

Removed deprecated demos folder.

parent 3e8668ee
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 875 deletions
nifty/demos/__init__.py deleted 100644 → 0
View file @ 3e8668ee
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2014 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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/>.
from demos_core import get_demo_dir
nifty/demos/demo_excaliwir.py deleted 100644 → 0
View file @ 3e8668ee
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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/>.
"""
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / demo
.. /______/
NIFTY demo for (extended) critical Wiener filtering of Gaussian random signals.
"""
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 *
##=============================================================================
class problem(object):
def __init__(self, x_space, s2n=6, **kwargs):
"""
Sets up a Wiener filter problem.
Parameters
----------
x_space : space
Position space the signal lives in.
s2n : float, *optional*
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.power_indices.set_default()
#self.k.set_power_indices(**kwargs)
## set some power spectrum
self.power = (lambda k: 42 / (k + 1) ** 2)
## define signal covariance
self.S = power_operator(self.k, spec=self.power, bare=True)
## define projector to spectral bands
self.Sk = self.S.get_projection_operator()
## generate signal
self.s = self.S.get_random_field(domain=self.z)
## define response
self.R = response_operator(self.z, sigma=0.0, mask=1.0)
## get data space
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), bare=True)
## define (plain) projector
self.Nj = projection_operator(d_space)
## generate noise
n = self.N.get_random_field(domain=d_space)
## compute data
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))
## define information propagator
self.D = propagator_operator(S=self.S,
N=self.N,
R=self.R)
## reserve map
self.m = None
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def solve(self, newspec=None):
"""
Solves the Wiener filter problem for a given power spectrum
reconstructing a signal estimate.
Parameters
----------
newspace : {scalar, list, array, Field, function}, *optional*
Assumed power spectrum (default: k ** -2).
"""
## set (given) power spectrum
if(newspec is None):
newspec = np.r_[1, 1 / self.k.power_indices["kindex"][1:] ** 2] ## Laplacian
elif(newspec is False):
newspec = self.power ## assumed to be known
self.S.set_power(newspec, bare=True)
## reconstruct map
self.m = self.D(self.j, W=self.S, tol=1E-3, note=False)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def solve_critical(self, newspec=None, q=0, alpha=1, delta=1, epsilon=0):
"""
Solves the (generalised) Wiener filter problem
reconstructing a signal estimate and a power spectrum.
Parameters
----------
newspace : {scalar, list, array, Field, function}, *optional*
Initial power spectrum (default: k ** -2).
q : {scalar, list, array}, *optional*
Spectral scale parameter of the assumed inverse-Gamme prior
(default: 0).
alpha : {scalar, list, array}, *optional*
Spectral shape parameter of the assumed inverse-Gamme prior
(default: 1).
delta : float, *optional*
First filter perception parameter (default: 1).
epsilon : float, *optional*
Second filter perception parameter (default: 0).
See Also
--------
infer_power
"""
## set (initial) power spectrum
if(newspec is None):
newspec = np.r_[1, 1 / self.k.power_indices["kindex"][1:] ** 2] ## Laplacian
elif(newspec is False):
newspec = self.power ## assumed to be known
self.S.set_power(newspec, bare=True)
## 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=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)
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)
#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.
"""
## plot signal
self.s.plot(title="signal")
## plot data
try:
d_ = Field(self.z, val=self.d.val, target=self.k)
d_.plot(title="data", vmin=self.s.min(), vmax=self.s.max())
except:
pass
## plot map
if(self.m is None):
self.s.plot(power=True, mono=False, other=self.power)
else:
self.m.plot(title="reconstructed map", vmin=self.s.min(), vmax=self.s.max())
self.m.plot(power=True, mono=False, other=(self.power, self.S.get_power()))
##=============================================================================
##-----------------------------------------------------------------------------
#
if(__name__=="__main__"):
x = RGSpace((1280), zerocenter=True)
p = problem(x, log = False)
## 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
##-----------------------------------------------------------------------------
nifty/demos/demo_faraday.py deleted 100644 → 0
View file @ 3e8668ee
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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/>.
"""
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / demo
.. /______/
NIFTY demo applying transformations to the "Faraday Map" [#]_.
References
----------
.. [#] N. Opermann et. al., "An improved map of the Galactic Faraday sky",
Astronomy & Astrophysics, Volume 542, id.A93, 06/2012;
`arXiv:1111.6186 <http://www.arxiv.org/abs/1111.6186>`_
"""
from __future__ import division
from nifty import *
##-----------------------------------------------------------------------------
# (global) Faraday map
m = field(HPSpace(128), val=np.load(os.path.join(get_demo_dir(),
"demo_faraday_map.npy")))
##-----------------------------------------------------------------------------
##=============================================================================
def run(projection=False, power=False):
"""
Runs the demo.
Parameters
----------
projection : bool, *optional*
Whether to additionaly show projections or not (default: False).
power : bool, *optional*
Whether to additionaly show power spectra or not (default: False).
"""
nicely = {"vmin":-4, "vmax":4, "cmap":ncmap.fm()}
# (a) representation on HEALPix grid
m0 = m
m0.plot(title=r"$m$ on a HEALPix grid", **nicely)
nicely.update({"cmap":ncmap.fm()}) # healpy bug workaround
# (b) representation in spherical harmonics
k_space = m0.target # == lm_space(383, mmax=383)
m1 = m0.transform(k_space) # == m.transform()
# m1.plot(title=r"$m$ in spherical harmonics")
if(power):
m1.plot(title=r"angular power spectrum of $m$", vmin=1E-2, vmax=1E+1, mono=False)
if(projection):
# define projection operator
Sk = projection_operator(m1.domain)
# project quadrupole
m2 = Sk(m0, band=2)
m2.plot(title=r"angular quadrupole of $m$ on a HEALPix grid", **nicely)
# (c) representation on Gauss-Legendre grid
y_space = m.target.get_codomain(coname="gl") # == gl_space(384, nlon=767)
m3 = m1.transform(y_space) # == m0.transform().transform(y_space)
m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", **nicely)
if(projection):
m4 = Sk(m3, band=2)
m4.plot(title=r"angular quadrupole of $m$ on a Gauss-Legendre grid", **nicely)
# (d) representation on regular grid
y_space = GLSpace(384, nlon=768) # auxiliary gl_space
z_space = RGSpace([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.paradict['nlon']//2, axis=1)) # rearrange value array
m5.plot(title=r"$m$ on a regular 2D grid", **nicely)
if(power):
m5.target.set_power_indices(log=False)
m5.plot(power=True, title=r"Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False)
if(projection):
m5.target.set_power_indices(log=False)
# define projection operator
Sk = projection_operator(m5.target)
# project quadrupole
m6 = Sk(m5, band=2)
m6.plot(title=r"Fourier quadrupole of $m$ on a regular 2D grid", **nicely)
##=============================================================================
##-----------------------------------------------------------------------------
if(__name__=="__main__"):
# pl.close("all")
# run demo
run(projection=False, power=False)
# define projection operator
Sk = projection_operator(m.target)
##-----------------------------------------------------------------------------
nifty/demos/demo_tomography.py deleted 100644 → 0
View file @ 3e8668ee
# -*- coding: utf-8 -*-
from nifty import *
if __name__ == "__main__":
shape = (256, 256)
x_space = RGSpace(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
nifty/demos/demo_wf1.py deleted 100644 → 0
View file @ 3e8668ee
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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/>.
"""
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / demo
.. /______/
NIFTY demo applying a Wiener filter using conjugate gradient.
"""
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__":
# some signal space; e.g., a two-dimensional regular grid
shape = [1024, 1024]
x_space = rg_space(shape)
#y_space = point_space(1280*1280)
#x_space = HPSpace(32)
#x_space = GLSpace(800)
k_space = x_space.get_codomain() # get conjugate space
# some power spectrum
power = (lambda k: 42 / (k + 1) ** 4)
S = power_operator(k_space, codomain=x_space, spec=power) # define signal covariance
s = S.get_random_field(domain=x_space) # generate signal
#my_mask = x_space.cast(1)
#stretch = 0.6
#my_mask[shape[0]/2*stretch:shape[0]/2/stretch, shape[1]/2*stretch:shape[1]/2/stretch] = 0
my_mask = 1
R = response_operator(x_space, sigma=0.01, mask=my_mask, assign=None) # define response
R = response_operator(x_space, assign=None) #
#R = identity_operator(x_space)
d_space = R.target # get data space
# some noise variance; e.g., signal-to-noise ratio of 1
N = diagonal_operator(d_space, diag=s.var(), bare=True) # define noise covariance
n = N.get_random_field(domain=d_space) # generate noise
d = R(s) + n # compute data
j = R.adjoint_times(N.inverse_times(d)) # define information source
D = propagator_operator(S=S, N=N, R=R) # define information propagator
#m = D(j, W=S, tol=1E-8, limii=100, note=True)
#m = D(j, tol=1E-8, limii=20, note=True, force=True)
ident = identity(x_space)
#xi = Field(x_space, random='gau', target=k_space)
m = D(j, W=S, tol=1E-8, limii=100, note=True)
#temp_result = (D.inverse_times(m)-xi)