Commit 39cdac31 by Martin Reinecke

### merge from master

parents b7fb0b3f 7c50e35f
 from nifty import * from nifty.library.wiener_filter import WienerFilterEnergy from nifty.library.critical_filter import CriticalPowerEnergy ... ... @@ -113,7 +112,7 @@ if __name__ == "__main__": print (x, iteration) minimizer1 = RelaxedNewton(convergence_tolerance=10e-2, minimizer1 = RelaxedNewton(convergence_tolerance=1e-2, convergence_level=2, iteration_limit=3, callback=convergence_measure) ... ... @@ -124,7 +123,7 @@ if __name__ == "__main__": max_history_length=3) # Setting starting position flat_power = Field(p_space,val=10e-8) flat_power = Field(p_space,val=1e-8) m0 = flat_power.power_synthesize(real_signal=True) t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2)) ... ... @@ -136,7 +135,7 @@ if __name__ == "__main__": distribution_strategy=distribution_strategy) # Initializing the nonlinear Wiener Filter energy map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=inverter) map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0) # Solving the Wiener Filter analytically D0 = map_energy.curvature m0 = D0.inverse_times(j) ... ...
 from nifty import * import plotly.offline as pl import plotly.graph_objs as go from nifty.library.wiener_filter import * from mpi4py import MPI comm = MPI.COMM_WORLD ... ... @@ -103,7 +103,7 @@ if __name__ == "__main__": # inverter = ConjugateGradient(convergence_level=3, convergence_tolerance=10e-5, convergence_tolerance=1e-5, preconditioner=None) # Setting starting position m0 = Field(h_space, val=.0) ... ...
 import numpy as np from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\ SmoothingOperator, DiagonalOperator, create_power_operator from nifty.library import WienerFilterCurvature from nifty import * #import plotly.offline as pl #import plotly.graph_objs as go ... ... @@ -10,36 +14,37 @@ 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 distribution_strategy = 'fftw' # 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) # 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) # variance of field in position space sqrt(<|s_x|^2>) (in unit of s) field_variance = 2. #smoothing length of response (in same unit as L) # smoothing length of response (in same unit as L) response_sigma = 0.1 #defining resolution (pixels per dimension) # defining resolution (pixels per dimension) N_pixels = 512 #Setting up derived constants # 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 # 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 pixel_length = L/N_pixels # Setting up the geometry s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width) fft = FFTOperator(s_space) s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length) fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex) h_space = fft.target[0] inverse_fft = FFTOperator(h_space, target=s_space, domain_dtype=np.complex, target_dtype=np.float) p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy) # Creating the mock data S = create_power_operator(h_space, power_spectrum=pow_spec, ... ... @@ -51,6 +56,7 @@ if __name__ == "__main__": ss = fft.inverse_times(sh) R = SmoothingOperator(s_space, sigma=response_sigma) R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0]) signal_to_noise = 1 N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True) ... ... @@ -63,7 +69,9 @@ if __name__ == "__main__": # Wiener filter j = R.adjoint_times(N.inverse_times(d)) D = PropagatorOperator(S=S, N=N, R=R) j = R_harmonic.adjoint_times(N.inverse_times(d)) wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic) m = wiener_curvature.inverse_times(j) m_s = inverse_fft(m) m = D(j)
 ... ... @@ -66,7 +66,7 @@ variable_default_field_dtype = keepers.Variable( variable_default_distribution_strategy = keepers.Variable( 'default_distribution_strategy', ['not', 'fftw', 'equal'], lambda z: (('pyfftw' in dependency_injector) lambda z: (('fftw' in dependency_injector) if z == 'fftw' else True), genus='str') ... ...
 ... ... @@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object): @staticmethod def _hermitian_decomposition(domain, val, spaces, domain_axes, preserve_gaussian_variance=False): # 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 spaces[1:]: (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() full = (hh + ha + ah + aa) h = (full + c)/2. a = (full - c)/2. flipped_val = val for space in spaces: flipped_val = domain[space].hermitianize_inverter( x=flipped_val, axes=domain_axes[space]) flipped_val = flipped_val.conjugate() h = (val + flipped_val)/2. a = val - h # correct variance if preserve_gaussian_variance: ... ...
 ... ... @@ -117,7 +117,6 @@ class FFTOperator(LinearOperator): super(FFTOperator, self).__init__(default_spaces) # Initialize domain and target self._domain = self._parse_domain(domain) if len(self.domain) != 1: raise ValueError("TransformationOperator accepts only exactly one " ... ...
 ... ... @@ -8,12 +8,21 @@ from .smoothing_operator import SmoothingOperator class FFTSmoothingOperator(SmoothingOperator): def _smooth(self, x, spaces, inverse): Transformator = FFTOperator(x.domain[spaces[0]]) def __init__(self, domain, sigma, log_distances=False, default_spaces=None): super(FFTSmoothingOperator, self).__init__( domain=domain, sigma=sigma, log_distances=log_distances, default_spaces=default_spaces) self._transformator_cache = {} def _smooth(self, x, spaces, inverse): # transform to the (global-)default codomain and perform all remaining # steps therein transformed_x = Transformator(x, spaces=spaces) transformator = self._get_transformator(x.dtype) transformed_x = transformator(x, spaces=spaces) codomain = transformed_x.domain[spaces[0]] coaxes = transformed_x.domain_axes[spaces[0]] ... ... @@ -51,9 +60,18 @@ class FFTSmoothingOperator(SmoothingOperator): transformed_x.val.set_local_data(local_transformed_x, copy=False) smoothed_x = Transformator.adjoint_times(transformed_x, spaces=spaces) smoothed_x = transformator.adjoint_times(transformed_x, spaces=spaces) result = x.copy_empty() result.set_val(smoothed_x, copy=False) return result def _get_transformator(self, dtype): if dtype not in self._transformator_cache: self._transformator_cache[dtype] = FFTOperator( self.domain, domain_dtype=dtype, target_dtype=np.complex) return self._transformator_cache[dtype]
 ... ... @@ -89,22 +89,6 @@ class LMSpace(Space): super(LMSpace, self).__init__() self._lmax = self._parse_lmax(lmax) def hermitian_decomposition(self, x, axes=None): if issubclass(x.dtype.type, np.complexfloating): hermitian_part = x.copy_empty() anti_hermitian_part = x.copy_empty() hermitian_part[:] = x.real anti_hermitian_part[:] = x.imag * 1j else: hermitian_part = x.copy() anti_hermitian_part = x.copy_empty() anti_hermitian_part[:] = 0 return (hermitian_part, anti_hermitian_part) def hermitian_fixed_points(self): return None # ---Mandatory properties and methods--- def __repr__(self): ... ...
 ... ... @@ -300,6 +300,9 @@ class PowerIndices(object): (temp_uniqued_pindex, local_temp_pundex) = np.unique( local_pindex, return_index=True) # cure a bug in numpy # https://github.com/numpy/numpy/issues/9405 local_temp_pundex = np.asarray(local_temp_pundex, dtype=np.int) # Shift the local pundices by the nodes' local_dim_offset local_temp_pundex += global_pindex.distributor.local_dim_offset ... ...
 ... ... @@ -91,26 +91,6 @@ class RGSpace(Space): self._shape = self._parse_shape(shape) self._distances = self._parse_distances(distances) def hermitian_decomposition(self, x, axes=None, preserve_gaussian_variance=False): # check axes if axes is None: axes = range(len(self.shape)) assert len(x.shape) >= len(self.shape), "shapes mismatch" assert len(axes) == len(self.shape), "axes mismatch" # compute the hermitian part flipped_x = self._hermitianize_inverter(x, axes=axes) flipped_x = flipped_x.conjugate() # average x and flipped_x. hermitian_part = x + flipped_x hermitian_part /= 2. # use subtraction since it is faster than flipping another time anti_hermitian_part = (x-hermitian_part) return (hermitian_part, anti_hermitian_part) def hermitian_fixed_points(self): dimensions = len(self.shape) mid_index = np.array(self.shape)//2 ... ... @@ -124,7 +104,7 @@ class RGSpace(Space): fixed_points += [tuple(index * mid_index)] return fixed_points def _hermitianize_inverter(self, x, axes): def hermitianize_inverter(self, x, axes): # calculate the number of dimensions the input array has dimensions = len(x.shape) # prepare the slicing object which will be used for mirroring ... ...
 ... ... @@ -147,35 +147,34 @@ class Space(DomainObject): raise NotImplementedError( "There is no generic co-smoothing kernel for Space base class.") def hermitian_decomposition(self, x, axes, preserve_gaussian_variance=False): """ Decomposes x into its hermitian and anti-hermitian constituents. def hermitian_fixed_points(self): """ Returns the array points which remain invariant under the action of `hermitianize_inverter` This method decomposes a field's array x into its hermitian and antihermitian part, which corresponds to real and imaginary part in a corresponding harmonic partner space. This is an internal function that is mainly used for power-synthesizing and -analyzing Fields. Returns ------- list of index-tuples The list contains the index-coordinates of the invariant points. """ return None def hermitianize_inverter(self, x, axes): """ Inverts/flips x in the context of Hermitian decomposition. This method is mainly used for power-synthesizing and -analyzing Fields. Parameters ---------- x : distributed_data_object The field's val object. axes : tuple of ints Specifies the axes of x which correspond to this space. preserve_gaussian_variance : bool *optional* If the hermitian decomposition is done via computing the half sums and differences of `x` and mirrored `x`, all points except the fixed points lose half of their variance. If `x` is complex also they lose half of their variance since the real(/imaginary) part gets lost. Returns ------- (distributed_data_object, distributed_data_object) A tuple of two distributed_data_objects, the first being the hermitian and the second the anti-hermitian part of x. distributed_data_object The Hermitian-flipped of x. """ raise NotImplementedError return x
 ... ... @@ -20,7 +20,7 @@ import unittest import numpy as np from numpy.testing import assert_, assert_equal, assert_raises,\ assert_almost_equal assert_almost_equal, assert_array_almost_equal from d2o import distributed_data_object from nifty import LMSpace from test.common import expand ... ... @@ -108,15 +108,12 @@ class LMSpaceFunctionalityTests(unittest.TestCase): for key, value in expected.iteritems(): assert_equal(getattr(l, key), value) @expand(get_hermitian_configs()) def test_hermitian_decomposition(self, x, real, imag): def test_hermitianize_inverter(self): l = LMSpace(5) assert_almost_equal( l.hermitian_decomposition(distributed_data_object(x))[0], real) assert_almost_equal( l.hermitian_decomposition(distributed_data_object(x))[1], imag*1j) v = distributed_data_object(global_shape=l.shape, dtype=np.complex128) v[:] = np.random.random(l.shape) + 1j*np.random.random(l.shape) inverted = l.hermitianize_inverter(v, axes=(0,)) assert_array_almost_equal(inverted.get_full_data(), v.get_full_data()) @expand(get_weight_configs()) def test_weight(self, x, power, axes, inplace, expected): ... ...
 ... ... @@ -143,7 +143,8 @@ class PowerSpaceConsistencyCheck(unittest.TestCase): distribution_strategy=distribution_strategy, logarithmic=logarithmic, nbin=nbin, binbounds=binbounds) assert_equal(p.pindex.flatten()[p.pundex], np.arange(p.dim), assert_equal(p.pindex.flatten().get_full_data()[p.pundex], np.arange(p.dim), err_msg='pundex is not right-inverse of pindex!') @expand(CONSISTENCY_CONFIGS) ... ... @@ -152,7 +153,6 @@ class PowerSpaceConsistencyCheck(unittest.TestCase): logarithmic): if distribution_strategy == "fftw": if not hasattr(gdi.get('fftw'), 'FFTW_MPI'): print (gdi.get('fftw'), "blub \n\n\n") raise SkipTest p = PowerSpace(harmonic_partner=harmonic_partner, distribution_strategy=distribution_strategy, ... ...
 ... ... @@ -21,6 +21,8 @@ from __future__ import division import unittest import numpy as np from d2o import distributed_data_object from numpy.testing import assert_, assert_equal, assert_almost_equal from nifty import RGSpace from test.common import expand ... ... @@ -127,49 +129,21 @@ class RGSpaceFunctionalityTests(unittest.TestCase): @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), (4, 6, 8), (17, 5, 3)])) def test_hermitian_decomposition(self, shape): def test_hermitianize_inverter(self, shape): r = RGSpace(shape, harmonic=True) v = np.empty(shape, dtype=np.complex128) v.real = np.random.random(shape) v.imag = np.random.random(shape) h, a = r.hermitian_decomposition(v) # make sure that data == h + a # NOTE: this is only correct for preserve_gaussian_variance==False, # but I consider this an intrinsic property of a hermitian # decomposition. assert_almost_equal(v, h+a) print (h, a) # test hermitianity of h it = np.nditer(h, flags=['multi_index']) while not it.finished: i1 = it.multi_index i2 = [] for i in range(len(i1)): i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0) i2 = tuple(i2) assert_almost_equal(h[i1], np.conj(h[i2])) assert_almost_equal(a[i1], -np.conj(a[i2])) it.iternext() v = distributed_data_object(global_shape=shape, dtype=np.complex128) v[:] = np.random.random(shape) + 1j*np.random.random(shape) inverted = r.hermitianize_inverter(v, axes=range(len(shape))) @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), (4, 6, 8), (17, 5, 3)])) def test_hermitian_decomposition2(self, shape): r = RGSpace(shape, harmonic=True) v = np.random.random(shape) h, a = r.hermitian_decomposition(v) # make sure that data == h + a assert_almost_equal(v, h+a) # test hermitianity of h it = np.nditer(h, flags=['multi_index']) # test hermitian flipping of `inverted` it = np.nditer(v, flags=['multi_index']) while not it.finished: i1 = it.multi_index i2 = [] for i in range(len(i1)): i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0) i2.append(v.shape[i]-i1[i] if i1[i] > 0 else 0) i2 = tuple(i2) assert_almost_equal(h[i1], np.conj(h[i2])) assert_almost_equal(a[i1], -np.conj(a[i2])) assert_almost_equal(inverted[i1], v[i2]) it.iternext() @expand(get_distance_array_configs()) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!