From ef955b6908ce6b4fff550dc3303b51a2eeaa5cee Mon Sep 17 00:00:00 2001 From: Ultima <theo.steininger@ultimanet.de> Date: Thu, 12 Nov 2015 09:22:55 +0100 Subject: [PATCH] Some Bugfixes. propagator_operator still not working to full extend. --- demos/demo_excaliwir.py | 15 ++- demos/demo_tomography.py | 33 +++--- demos/demo_wf1.py | 51 +++++--- nifty_core.py | 8 +- nifty_mpi_data.py | 1 + nifty_power_indices.py | 9 +- nifty_utilities.py | 82 ++++++++----- operators/nifty_operators.py | 17 +-- operators/nifty_probing.py | 220 +++++++++++++---------------------- rg/nifty_rg.py | 53 ++++----- 10 files changed, 246 insertions(+), 243 deletions(-) diff --git a/demos/demo_excaliwir.py b/demos/demo_excaliwir.py index c3281052a..8d25f1912 100644 --- a/demos/demo_excaliwir.py +++ b/demos/demo_excaliwir.py @@ -46,7 +46,7 @@ from nifty import * class problem(object): - def __init__(self, x_space, s2n=0.5, **kwargs): + def __init__(self, x_space, s2n=12, **kwargs): """ Sets up a Wiener filter problem. @@ -67,7 +67,7 @@ class problem(object): #self.k.set_power_indices(**kwargs) ## set some power spectrum - self.power = (lambda k: 42 / (k + 1) ** 5) + self.power = (lambda k: 42 / (k + 1) ** 3) ## define signal covariance self.S = power_operator(self.k, spec=self.power, bare=True) @@ -82,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 @@ -95,7 +96,9 @@ class problem(object): #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 @@ -186,7 +189,7 @@ class problem(object): print ('power', power) print 'Calculated power' - power = np.clip(power, 0.00000001, np.max(power)) + #power = np.clip(power, 0.00000001, np.max(power)) self.store += [{'tr_B1': tr_B1, 'tr_B2': tr_B2, 'num': numerator, @@ -253,7 +256,7 @@ class problem(object): ##----------------------------------------------------------------------------- # if(__name__=="__main__"): - x = rg_space((128,128), zerocenter=True) + x = rg_space((128), zerocenter=True) p = problem(x, log = False) about.warnings.off() ## pl.close("all") diff --git a/demos/demo_tomography.py b/demos/demo_tomography.py index 8ef3cf593..f186f3610 100644 --- a/demos/demo_tomography.py +++ b/demos/demo_tomography.py @@ -8,12 +8,21 @@ shape = (256, 256) x_space = rg_space(shape) k_space = x_space.get_codomain() -power = lambda k: 42/((1+k*shape[0])**2) +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 = [] @@ -29,9 +38,9 @@ def make_los(n=10, angle=0, d=1): 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 = 256 -m = 256 +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) @@ -39,27 +48,19 @@ for alpha in [np.pi/n*j for j in xrange(n)]: starts = list(temp_coords[0].T) ends = list(temp_coords[1].T) -#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) - -R = los_response(x_space, starts=starts, ends=ends, sigmas_up=0.1, sigmas_low=0.1) +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(), bare=True) +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) -m = D(j, W=S, tol=1E-14, limii=50, note=True) s.plot(title="signal", save='1_plot_s.png') s.plot(save='plot_s_power.png', power=True, other=power) diff --git a/demos/demo_wf1.py b/demos/demo_wf1.py index 8a01d0072..5543a479f 100644 --- a/demos/demo_wf1.py +++ b/demos/demo_wf1.py @@ -44,25 +44,28 @@ from nifty import * # version 0.8.0 about.warnings.off() # some signal space; e.g., a two-dimensional regular grid -x_space = rg_space([128, 128]) # define signal space +shape = [1024,] +x_space = rg_space(shape) #y_space = point_space(1280*1280) - #x_space = hp_space(32) - #x_space = gl_space(800) k_space = x_space.get_codomain() # get conjugate space # some power spectrum -power = (lambda k: 42 / (k + 1) ** 3) +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) -#my_mask[400:900,400:900] = 0 +#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 @@ -79,14 +82,34 @@ D = propagator_operator(S=S, N=N, R=R) # define inform #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(xi, W=ident, tol=1E-8, limii=10, note=True, force=True) -temp_result = (D.inverse_times(m)-xi) -print (temp_result.dot(temp_result)) -print (temp_result.val) +#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) + -s.plot(title="signal", save = 'plot_s.png') # plot signal -d_ = field(x_space, val=d.val, target=k_space) -d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png') # plot data -m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png') # plot map +#n_power = x_space.enforce_power(s.var()/np.prod(shape)) +#s_power = S.get_power() + +#s.plot(title="signal", save = 'plot_s.png') +#s.plot(title="signal power", power=True, other=power, +# mono=False, save = 'power_plot_s.png', nbin=1000, log=True, +# vmax = 100, vmin=10e-7) + +#d_ = field(x_space, val=d.val, target=k_space) +#d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png') + + +#n_ = field(x_space, val=n.val, target=k_space) +#n_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_n.png') + + +# +#m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png') +#m.plot(title="reconstructed power", power=True, other=(n_power, s_power), +# save = 'power_plot_m.png', vmin=0.001, vmax=10, mono=False) +# # diff --git a/nifty_core.py b/nifty_core.py index e5e538a2f..d64305f09 100644 --- a/nifty_core.py +++ b/nifty_core.py @@ -1921,7 +1921,8 @@ class field(object): """ - def __init__(self, domain, val=None, codomain=None, ishape=None, **kwargs): + def __init__(self, domain=None, val=None, codomain=None, ishape=None, + **kwargs): """ Sets the attributes for a field class instance. @@ -1944,6 +1945,9 @@ class field(object): Nothing """ + # If the given val was a field, try to cast it accordingly to the given + # domain and codomain, etc... + # check domain if not isinstance(domain, space): raise TypeError(about._errors.cstring("ERROR: invalid input.")) @@ -2359,7 +2363,7 @@ class field(object): # whether the domain matches exactly or not: # extract the data from x and try to dot with this - return self.dot(x=x.get_val(), axis=axis) + return self.dot(x=x.get_val(), axis=axis, bare=bare) # Case 3: x is something else else: diff --git a/nifty_mpi_data.py b/nifty_mpi_data.py index 559096421..8790d30cd 100644 --- a/nifty_mpi_data.py +++ b/nifty_mpi_data.py @@ -631,6 +631,7 @@ class distributed_data_object(object): temp_d2o = self.copy_empty() temp_data = np.conj(self.get_local_data()) temp_d2o.set_local_data(temp_data) + temp_d2o.hermitian = self.hermitian return temp_d2o def conj(self): diff --git a/nifty_power_indices.py b/nifty_power_indices.py index c7b2b5b5c..c5cc5f76e 100644 --- a/nifty_power_indices.py +++ b/nifty_power_indices.py @@ -53,7 +53,9 @@ class power_indices(object): # Initialize the dictonary which stores all individual index-dicts self.global_dict = {} # Set self.default_parameters - self.set_default(log=log, nbin=nbin, binbounds=binbounds) + self.set_default(config_dict={'log': log, + 'nbin': nbin, + 'binbounds': binbounds}) # Redirect the direct calls approaching a power_index instance to the # default_indices dict @@ -217,7 +219,10 @@ class power_indices(object): # indices, bin them, compute the pundex and then return everything. else: # Get the unbinned indices - temp_unbinned_indices = self.get_index_dict(store=False) + temp_unbinned_indices = self.get_index_dict(nbin=None, + binbounds=None, + log=False, + store=False) # Bin them (temp_pindex, temp_kindex, temp_rho, temp_pundex) = \ self._bin_power_indices( diff --git a/nifty_utilities.py b/nifty_utilities.py index b40b063fb..7519ecc9e 100644 --- a/nifty_utilities.py +++ b/nifty_utilities.py @@ -2,23 +2,24 @@ import numpy as np -def hermitianize(x): - ## make the point inversions + +def hermitianize_gaussian(x): + # make the point inversions flipped_x = _hermitianize_inverter(x) flipped_x = flipped_x.conjugate() - ## check if x was already hermitian + # check if x was already hermitian if (x == flipped_x).all(): return x - ## average x and flipped_x. - ## Correct the variance by multiplying sqrt(0.5) + # average x and flipped_x. + # Correct the variance by multiplying sqrt(0.5) x = (x + flipped_x) * np.sqrt(0.5) - ## The fixed points of the point inversion must not be avaraged. - ## Hence one must multiply them again with sqrt(0.5) - ## -> Get the middle index of the array + # The fixed points of the point inversion must not be avaraged. + # Hence one must multiply them again with sqrt(0.5) + # -> Get the middle index of the array mid_index = np.array(x.shape, dtype=np.int)//2 dimensions = mid_index.size - ## Use ndindex to iterate over all combinations of zeros and the - ## mid_index in order to correct all fixed points. + # Use ndindex to iterate over all combinations of zeros and the + # mid_index in order to correct all fixed points. for i in np.ndindex((2,)*dimensions): temp_index = tuple(i*mid_index) x[temp_index] *= np.sqrt(0.5) @@ -30,18 +31,35 @@ def hermitianize(x): return x +def hermitianize(x): + # make the point inversions + flipped_x = _hermitianize_inverter(x) + flipped_x = flipped_x.conjugate() + # check if x was already hermitian + if (x == flipped_x).all(): + return x + # average x and flipped_x. + # Correct the variance by multiplying sqrt(0.5) + x = (x + flipped_x) / 2. + try: + x.hermitian = True + except(AttributeError): + pass + + return x + def _hermitianize_inverter(x): - ## calculate the number of dimensions the input array has + # calculate the number of dimensions the input array has dimensions = len(x.shape) - ## prepare the slicing object which will be used for mirroring - slice_primitive = [slice(None),]*dimensions - ## copy the input data + # prepare the slicing object which will be used for mirroring + slice_primitive = [slice(None), ]*dimensions + # copy the input data y = x.copy() - ## flip in every direction + # flip in every direction for i in xrange(dimensions): slice_picker = slice_primitive[:] - slice_picker[i] = slice(1, None,None) + slice_picker[i] = slice(1, None, None) slice_inverter = slice_primitive[:] slice_inverter[i] = slice(None, 0, -1) @@ -54,7 +72,7 @@ def _hermitianize_inverter(x): def direct_dot(x, y): - ## the input could be fields. Try to extract the data + # the input could be fields. Try to extract the data try: x = x.get_val() except(AttributeError): @@ -63,7 +81,7 @@ def direct_dot(x, y): y = y.get_val() except(AttributeError): pass - ## try to make a direct vdot + # try to make a direct vdot try: return x.vdot(y) except(AttributeError): @@ -74,16 +92,16 @@ def direct_dot(x, y): except(AttributeError): pass - ## fallback to numpy + # fallback to numpy return np.vdot(x, y) def convert_nested_list_to_object_array(x): - ## if x is a nested_list full of ndarrays all having the same size, - ## np.shape returns the shape of the ndarrays, too, i.e. too many - ## dimensions + # if x is a nested_list full of ndarrays all having the same size, + # np.shape returns the shape of the ndarrays, too, i.e. too many + # dimensions possible_shape = np.shape(x) - ## Check if possible_shape goes too deep. + # Check if possible_shape goes too deep. dimension_counter = 0 current_extract = x for i in xrange(len(possible_shape)): @@ -93,11 +111,11 @@ def convert_nested_list_to_object_array(x): current_extract = current_extract[0] dimension_counter += 1 real_shape = possible_shape[:dimension_counter] - ## if the numpy array was not encapsulated at all, return x directly + # if the numpy array was not encapsulated at all, return x directly if real_shape == (): return x - ## Prepare the carrier-object - carrier = np.empty(real_shape, dtype = np.object) + # Prepare the carrier-object + carrier = np.empty(real_shape, dtype=np.object) for i in xrange(np.prod(real_shape)): ii = np.unravel_index(i, real_shape) try: @@ -121,10 +139,10 @@ def field_map(ishape, function, *args): result[ii] = function() return result else: - ## define a helper function in order to clip the get-indices - ## to be suitable for the foreign arrays in args. - ## This allows you to do operations, like adding to fields - ## with ishape (3,4,3) and (3,4,1) + # define a helper function in order to clip the get-indices + # to be suitable for the foreign arrays in args. + # This allows you to do operations, like adding to fields + # with ishape (3,4,3) and (3,4,1) def get_clipped(w, ind): w_shape = np.array(np.shape(w)) get_tuple = tuple(np.clip(ind, 0, w_shape-1)) @@ -135,5 +153,5 @@ def field_map(ishape, function, *args): result[ii] = function(*map( lambda z: get_clipped(z, ii), args) ) - #result[ii] = function(*map(lambda z: z[ii], args)) - return result \ No newline at end of file + # result[ii] = function(*map(lambda z: z[ii], args)) + return result diff --git a/operators/nifty_operators.py b/operators/nifty_operators.py index af1fcab05..6b800b170 100644 --- a/operators/nifty_operators.py +++ b/operators/nifty_operators.py @@ -2292,7 +2292,7 @@ class projection_operator(operator): raise TypeError(about._errors.cstring("ERROR: Invalid bands.")) if bands_was_scalar: - new_field = x * (self.assign == bands[0]) + new_field = fx * (self.assign == bands[0]) else: # build up the projection results # prepare the projector-carrier @@ -2392,15 +2392,18 @@ class projection_operator(operator): vecvec = vecvec_operator(val=x) return self.pseudo_tr(x=vecvec, axis=axis, **kwargs) - # Case 2: x is not an operator - elif isinstance(x, operator) == False: + # Case 2: x is an operator + # -> take the diagonal + elif isinstance(x, operator): + working_field = x.diag(bare=False) + if self.domain != working_field.domain: + working_field = working_field.transform(codomain=self.domain) + + # Case 3: x is something else + else: raise TypeError(about._errors.cstring( "ERROR: x must be a field or an operator.")) - # Case 3: x is an operator - # -> take the diagonal - working_field = x.diag() - # Check for hidden degrees of freedom and compensate the trace # accordingly if self.domain.get_dim() != self.domain.get_dof(): diff --git a/operators/nifty_probing.py b/operators/nifty_probing.py index 8b75db42a..aee9d55e4 100644 --- a/operators/nifty_probing.py +++ b/operators/nifty_probing.py @@ -1,23 +1,23 @@ -## NIFTY (Numerical Information Field Theory) has been developed at the -## Max-Planck-Institute for Astrophysics. -## -## Copyright (C) 2015 Max-Planck-Society -## -## Author: Theo Steininger -## 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/>. +# NIFTY (Numerical Information Field Theory) has been developed at the +# Max-Planck-Institute for Astrophysics. +# +# Copyright (C) 2015 Max-Planck-Society +# +# Author: Theo Steininger +# 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 __future__ import division @@ -27,7 +27,7 @@ from nifty.nifty_core import space, \ from nifty.nifty_utilities import direct_dot -##============================================================================= + class prober(object): """ @@ -177,59 +177,61 @@ class prober(object): zeroth entry and the variance in the first entry. (default: False) """ - ## Case 1: no operator given. Check function and domain for general - ## sanity + # Case 1: no operator given. Check function and domain for general + # sanity if operator is None: - ## check whether the given function callable - if function is None or hasattr(function, "__call__") == False: + # check whether the given function callable + if function is None or not hasattr(function, "__call__"): raise ValueError(about._errors.cstring( - "ERROR: invalid input: No function given or not callable.")) - ## check given domain - if domain is None or isinstance(domain, space) == False: + "ERROR: invalid input: No function given or not callable.")) + # check given domain + if domain is None or not isinstance(domain, space): raise ValueError(about._errors.cstring( - "ERROR: invalid input: given domain is not a nifty space")) + "ERROR: invalid input: given domain is not a nifty space")) - ## Case 2: An operator is given. Take domain and function from that - ## if not given explicitly + # Case 2: An operator is given. Take domain and function from that + # if not given explicitly else: - ## Case 2.1 extract function - ## explicit function overrides operator function - if function is None or hasattr(function,"__call__") == False: + # Check 2.1 extract function + # explicit function overrides operator function + if function is None or not hasattr(function, "__call__"): try: function = operator.times except(AttributeError): raise ValueError(about._errors.cstring( - "ERROR: no explicit function given and given operator has no times method!")) - ## check whether the given function is correctly bound to the - ## operator + "ERROR: no explicit function given and given " + + "operator has no times method!")) + + # Check 2.2 check whether the given function is correctly bound to + # the operator if operator != function.im_self: raise ValueError(about._errors.cstring( - "ERROR: the given function is not a bound function of the operator!")) - - - - ## Case 2.2 extract domain - if domain is None or isinstance(domain, space): - if (function in [operator.inverse_times, - operator.adjoint_times]): - try: - domain = operator.target - except(AttributeError): - raise ValueError(about._errors.cstring( - "ERROR: no explicit domain given and given operator has no target!")) - - else: - try: - domain = operator.domain - except(AttributeError): - raise ValueError(about._errors.cstring( - "ERROR: no explicit domain given and given operator has no domain!")) + "ERROR: the given function is not a bound function " + + "of the operator!")) + + # Check 2.3 extract domain + if domain is None or not isinstance(domain, space): + if (function in [operator.inverse_times, + operator.adjoint_times]): + try: + domain = operator.target + except(AttributeError): + raise ValueError(about._errors.cstring( + "ERROR: no explicit domain given and given " + + "operator has no target!")) + else: + try: + domain = operator.domain + except(AttributeError): + raise ValueError(about._errors.cstring( + "ERROR: no explicit domain given and given " + + "operator has no domain!")) self.function = function self.domain = domain - ## Check the given target + # Check the given codomain if codomain is None: codomain = self.domain.get_codomain() else: @@ -237,71 +239,16 @@ class prober(object): self.codomain = codomain - if(random not in ["pm1","gau"]): + if random not in ["pm1", "gau"]: raise ValueError(about._errors.cstring( - "ERROR: unsupported random key '"+str(random)+"'.")) + "ERROR: unsupported random key '" + str(random) + "'.")) self.random = random - ## Parse the remaining arguments + # Parse the remaining arguments self.nrun = int(nrun) self.varQ = bool(varQ) self.kwargs = kwargs - """ - from nifty_operators import operator - if(not isinstance(op,operator)): - raise TypeError(about._errors.cstring("ERROR: invalid input.")) - ## check whether callable - if(function is None)or(not hasattr(function,"__call__")): - function = op.times - elif(op==function): - function = op.times - ## check whether correctly bound - if(op!=function.im_self): - raise NameError(about._errors.cstring("ERROR: invalid input.")) - ## check given shape and domain - if(domain is None)or(not isinstance(domain,space)): - if(function in [op.inverse_times,op.adjoint_times]): - domain = op.target - else: - domain = op.domain - else: - if(function in [op.inverse_times,op.adjoint_times]): - op.target.check_codomain(domain) ## a bit pointless - if(target is None)or(not isinstance(target,space)): - target = op.target - else: - op.domain.check_codomain(domain) ## a bit pointless - if(target is None)or(not isinstance(target,space)): - target = op.domain - - self.function = function - self.domain = domain - - ## check codomain - if(target is None): - target = self.domain.get_codomain() - else: - self.domain.check_codomain(target) ## a bit pointless - self.target = target - - if(random not in ["pm1","gau"]): - raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(random)+"'.")) - self.random = random - - self.ncpu = int(max(1,ncpu)) - self.nrun = int(max(self.ncpu,nrun)) - if(nper is None): - self.nper = None - else: - self.nper = int(max(1,min(self.nrun//self.ncpu,nper))) - - self.var = bool(var) - - self.quargs = quargs - """ - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - def configure(self, **kwargs): """ changes the attributes of the instance @@ -321,10 +268,11 @@ class prober(object): whether the variance will be additionally returned (default: False) """ - if("random" in kwargs): - if kwargs.get("random") not in ["pm1","gau"]: + if "random" in kwargs: + if kwargs.get("random") not in ["pm1", "gau"]: raise ValueError(about._errors.cstring( - "ERROR: unsupported random key '"+str(kwargs.get("random"))+"'.")) + "ERROR: unsupported random key '" + + str(kwargs.get("random")) + "'.")) else: self.random = kwargs.get("random") @@ -334,8 +282,6 @@ class prober(object): if "varQ" in kwargs: self.varQ = bool(kwargs.get("varQ")) - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - def generate_probe(self): """ Generates a single probe @@ -348,10 +294,8 @@ class prober(object): """ return field(self.domain, - codomain = self.codomain, - random = self.random) - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + codomain=self.codomain, + random=self.random) def evaluate_probe(self, probe, idnum=0): """ @@ -375,7 +319,7 @@ class prober(object): return f - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + def finalize(self, sum_of_probes, sum_of_squares, num): """ @@ -401,7 +345,7 @@ class prober(object): (`final`,`var`). """ - ## Check the success and efficiency of the probing + # Check the success and efficiency of the probing if num < self.nrun: about.infos.cflush( " ( %u probe(s) failed, effectiveness == %.1f%% )\n"\ @@ -424,22 +368,22 @@ class prober(object): else: return sum_of_probes*(1./num) - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - def print_progress(self, num): ## > prints progress status by in upto 10 dots + + def print_progress(self, num): # > prints progress status by in upto 10 dots tenths = 1+(10*num//self.nrun) about.infos.cflush(("\b")*10+('.')*tenths+(' ')*(10-tenths)) """ - def _single_probing(self,zipped): ## > performs one probing operation - ## generate probe + def _single_probing(self,zipped): # > performs one probing operation + # generate probe np.random.seed(zipped[0]) probe = self.gen_probe() - ## do the actual probing + # do the actual probing return self.probing(zipped[1],probe) """ - def probe(self): ## > performs the probing operations one after another - ## initialize the variables + def probe(self): # > performs the probing operations one after another + # initialize the variables sum_of_probes = 0 sum_of_squares = 0 num = 0 @@ -456,7 +400,7 @@ class prober(object): num += 1 self.print_progress(num) about.infos.cflush(" done.") - ## evaluate + # evaluate return self.finalize(sum_of_probes, sum_of_squares, num) def __call__(self,loop=False,**kwargs): @@ -484,18 +428,18 @@ class prober(object): self.configure(**kwargs) return self.probe() - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + def __repr__(self): return "<nifty_core.probing>" -##============================================================================= + class _specialized_prober(object): def __init__(self, operator, domain=None, inverseQ=False, **kwargs): - ## remove a potentially supplied function keyword argument + # remove a potentially supplied function keyword argument try: kwargs.pop('function') except(KeyError): diff --git a/rg/nifty_rg.py b/rg/nifty_rg.py index d89c8cdd7..9920ebc04 100644 --- a/rg/nifty_rg.py +++ b/rg/nifty_rg.py @@ -267,7 +267,7 @@ class rg_space(point_space): dtype=dtype, **kwargs) if x is not None and hermitianize and \ - self.paradict['complexity'] == 1 and not casted_x.hermitian: + self.paradict['complexity'] == 1 and not casted_x.hermitian: about.warnings.cflush( "WARNING: Data gets hermitianized. This operation is " + "extremely expensive\n") @@ -593,13 +593,18 @@ class rg_space(point_space): sample[sample >= 0] = 1 sample[sample < 0] = -1 + try: + sample.hermitian = True + except(AttributeError): + pass + # Case 2: normal distribution with zero-mean and a given standard # deviation or variance elif arg['random'] == 'gau': sample = super(rg_space, self).get_random_values(**arg) if hermitianizeQ: - sample = utilities.hermitianize(sample) + sample = utilities.hermitianize_gaussian(sample) # Case 3: uniform distribution elif arg['random'] == "uni" and not hermitianizeQ: @@ -633,6 +638,11 @@ class rg_space(point_space): sample *= (vmax - vmin) / 2. sample += 1 / 2. * (vmax + vmin) + try: + sample.hermitian = True + except(AttributeError): + pass + elif(arg['random'] == "syn"): spec = arg['spec'] kpack = arg['kpack'] @@ -714,11 +724,8 @@ class rg_space(point_space): elif self.datamodel in RG_DISTRIBUTION_STRATEGIES: # extract the local data from kdict local_kdict = kdict.get_local_data() - print ('local_kdict', local_kdict) rescaler = np.sqrt( spec[np.searchsorted(kindex, local_kdict)]) - print ('rescaler', rescaler) - print ('sample', sample.distribution_strategy) sample.apply_scalar_function(lambda x: x * rescaler, inplace=True) else: @@ -758,9 +765,7 @@ class rg_space(point_space): kindex=kpack[1], spec=spec, codomain=self, - log=log, - nbin=nbin, - binbounds=binbounds) + **lnb_dict) # Perform a fourier transform sample = temp_codomain.calc_transform(sample, codomain=self) @@ -1053,7 +1058,8 @@ class rg_space(point_space): power_spectrum /= rho return power_spectrum - def get_plot(self,x,title="",vmin=None,vmax=None,power=None,unit="",norm=None,cmap=None,cbar=True,other=None,legend=False,mono=True,**kwargs): + def get_plot(self,x,title="",vmin=None,vmax=None,power=None,unit="", + norm=None,cmap=None,cbar=True,other=None,legend=False,mono=True,**kwargs): """ Creates a plot of field values according to the specifications given by the parameters. @@ -1121,10 +1127,7 @@ class rg_space(point_space): (default: 0). """ - try: - x = x.get_full_data() - except AttributeError: - pass + if(not pl.isinteractive())and(not bool(kwargs.get("save",False))): about.warnings.cprint("WARNING: interactive mode off.") @@ -1134,8 +1137,13 @@ class rg_space(point_space): if(power): x = self.calc_power(x,**kwargs) + try: + x = x.get_full_data() + except AttributeError: + pass - fig = pl.figure(num=None,figsize=(6.4,4.8),dpi=None,facecolor="none",edgecolor="none",frameon=False,FigureClass=pl.Figure) + fig = pl.figure(num=None,figsize=(6.4,4.8),dpi=None, + facecolor="none",edgecolor="none",frameon=False,FigureClass=pl.Figure) ax0 = fig.add_axes([0.12,0.12,0.82,0.76]) ## explicit kindex @@ -1148,19 +1156,8 @@ class rg_space(point_space): except: kindex_supply_space = self.get_codomain() - - try: - default_indices = self.power_indices.default_parameters - except AttributeError: - default_indices = self.get_codomain().power_indices.default_parameters - log = kwargs.get('log', default_indices['log']) - nbin = kwargs.get('nbin', default_indices['nbin']) - binbounds = kwargs.get('binbounds', default_indices['binbounds']) - xaxes = kindex_supply_space.power_indices.get_index_dict( - log=log, - nbin=nbin, - binbounds=binbounds)['kindex'] + **kwargs)['kindex'] # try: @@ -1209,6 +1206,10 @@ class rg_space(point_space): ax0.set_title(title) else: + try: + x = x.get_full_data() + except AttributeError: + pass if(naxes==1): fig = pl.figure(num=None,figsize=(6.4,4.8),dpi=None,facecolor="none",edgecolor="none",frameon=False,FigureClass=pl.Figure) ax0 = fig.add_axes([0.12,0.12,0.82,0.76]) -- GitLab