diff --git a/__init__.py b/__init__.py index 2037d14f7658448b7d945b14b5f2ab9f6bb9212b..c33b568ba67743b8468f0337f9a7915a62c1e3d7 100644 --- a/__init__.py +++ b/__init__.py @@ -34,8 +34,7 @@ from nifty_mpi_data import distributed_data_object from nifty_power import * from nifty_random import random from nifty_simple_math import * -from nifty_tools import conjugate_gradient,\ - steepest_descent + from nifty_paradict import space_paradict,\ point_space_paradict,\ nested_space_paradict diff --git a/demos/demo_wf1.py b/demos/demo_wf1.py index ac96c0b0b1028e277a89973a38dae8c3b45b73b0..15439882d5248aa9285c3d907449af6fbe578b05 100644 --- a/demos/demo_wf1.py +++ b/demos/demo_wf1.py @@ -63,8 +63,8 @@ D = propagator_operator(S=S, N=N, R=R) # define inform m = D(j, W=S, tol=1E-3, note=True) # reconstruct map -s.plot(title="signal") # plot signal +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()) # plot data -m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max()) # plot map +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 diff --git a/lm/nifty_lm.py b/lm/nifty_lm.py index cab39354b7411d319e48b63fbb3e86db9d18af29..788f41d85b23433cfa2d6fd42886c44e74f51360 100644 --- a/lm/nifty_lm.py +++ b/lm/nifty_lm.py @@ -1372,6 +1372,13 @@ class gl_space(point_space): else: return gl.weight(x,self.vol,p=np.float64(power),nlat=self.para[0],nlon=self.para[1],overwrite=False) + def get_weight(self, power = 1): + ## TODO: Check if this function is compatible to the rest of the nifty code + ## TODO: Can this be done more efficiently? + dummy = self.enforce_values(1) + weighted_dummy = self.calc_weight(dummy, power = power) + return weighted_dummy/dummy + ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def calc_transform(self,x,codomain=None,**kwargs): diff --git a/nifty_core.py b/nifty_core.py index 3a8c72506813b9b7162e68784c30cb838af417c4..70fb33d02026f308d207e37fc3c2bc9cc7e4c86e 100644 --- a/nifty_core.py +++ b/nifty_core.py @@ -682,6 +682,9 @@ class space(object): """ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'calc_weight'.")) + def get_weight(self, power=1): + raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'get_weight'.")) + ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def calc_dot(self,x,y): @@ -1608,8 +1611,11 @@ class point_space(space): """ x = self.enforce_shape(np.array(x,dtype=self.datatype)) ## weight - return x*self.vol**power + return x*self.get_weight(power = power) + #return x*self.vol**power + def get_weight(self, power = 1): + return self.vol**power ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def calc_dot(self, x, y): """ @@ -2291,7 +2297,11 @@ class nested_space(space): """ x = self.enforce_shape(np.array(x,dtype=self.datatype)) ## weight - return x*self.get_meta_volume(total=False)**power + return x*self.get_weight(power = power) + + def get_weight(self, power = 1): + return self.get_meta_volume(total=False)**power + ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -2669,7 +2679,7 @@ class field(object): if val == None: if kwargs == {}: - self.val = self.domain.cast(0) + self.val = self.domain.cast(0.) else: self.val = self.domain.get_random_values(codomain=self.target, **kwargs) @@ -3349,8 +3359,9 @@ class field(object): temp = self else: temp = self.copy_empty() - data_object = self.domain.apply_scalar_function(self.val,\ - function, inplace) + data_object = self.domain.apply_scalar_function(self.val, + function, + inplace) temp.set_val(data_object) return temp diff --git a/nifty_mpi_data.py b/nifty_mpi_data.py index 47fdd5ae8e6854dbe1a6d3db9e63e1adfa306b5c..36821e031e2a4582f7468f3b290f7741f18978fa 100644 --- a/nifty_mpi_data.py +++ b/nifty_mpi_data.py @@ -163,11 +163,11 @@ class distributed_data_object(object): **kwargs) return temp_d2o - def apply_scalar_function(self, function, inplace=False): + def apply_scalar_function(self, function, inplace=False, dtype=None): if inplace == True: temp = self else: - temp = self.copy_empty() + temp = self.copy_empty(dtype=dtype) try: temp.data[:] = function(self.data) @@ -260,34 +260,54 @@ class distributed_data_object(object): temp_d2o.set_local_data(data = self.get_local_data().__abs__()) return temp_d2o - def __builtin_helper__(self, operator, other): + def __builtin_helper__(self, operator, other, inplace=False): ## Case 1: other is not a scalar if not (np.isscalar(other) or np.shape(other) == (1,)): ## if self.shape != other.shape: ## raise AttributeError(about._errors.cstring( ## "ERROR: Shapes do not match!")) - + try: + hermitian_Q = other.hermitian + except(AttributeError): + hermitian_Q = False ## extract the local data from the 'other' object temp_data = self.distributor.extract_local_data(other) temp_data = operator(temp_data) - else: + ## Case 2: other is a real scalar -> preserve hermitianity + elif np.isreal(other) or (self.dtype not in (np.complex, np.complex128, + np.complex256)): + hermitian_Q = self.hermitian temp_data = operator(other) - + ## Case 3: other is complex + else: + hermitian_Q = False + temp_data = operator(other) ## write the new data into a new distributed_data_object - temp_d2o = self.copy_empty() + if inplace == True: + temp_d2o = self + else: + temp_d2o = self.copy_empty() temp_d2o.set_local_data(data=temp_data) + temp_d2o.hermitian = hermitian_Q return temp_d2o - + """ def __inplace_builtin_helper__(self, operator, other): + ## Case 1: other is not a scalar if not (np.isscalar(other) or np.shape(other) == (1,)): temp_data = self.distributor.extract_local_data(other) temp_data = operator(temp_data) - else: + ## Case 2: other is a real scalar -> preserve hermitianity + elif np.isreal(other): + hermitian_Q = self.hermitian temp_data = operator(other) + ## Case 3: other is complex + else: + temp_data = operator(other) self.set_local_data(data=temp_data) + self.hermitian = hermitian_Q return self - + """ def __add__(self, other): return self.__builtin_helper__(self.get_local_data().__add__, other) @@ -296,8 +316,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__radd__, other) def __iadd__(self, other): - return self.__inplace_builtin_helper__(self.get_local_data().__iadd__, - other) + return self.__builtin_helper__(self.get_local_data().__iadd__, + other, + inplace = True) def __sub__(self, other): return self.__builtin_helper__(self.get_local_data().__sub__, other) @@ -306,8 +327,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__rsub__, other) def __isub__(self, other): - return self.__inplace_builtin_helper__(self.get_local_data().__isub__, - other) + return self.__builtin_helper__(self.get_local_data().__isub__, + other, + inplace = True) def __div__(self, other): return self.__builtin_helper__(self.get_local_data().__div__, other) @@ -316,8 +338,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__rdiv__, other) def __idiv__(self, other): - return self.__inplace_builtin_helper__(self.get_local_data().__idiv__, - other) + return self.__builtin_helper__(self.get_local_data().__idiv__, + other, + inplace = True) def __floordiv__(self, other): return self.__builtin_helper__(self.get_local_data().__floordiv__, @@ -326,8 +349,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__rfloordiv__, other) def __ifloordiv__(self, other): - return self.__inplace_builtin_helper__( - self.get_local_data().__ifloordiv__, other) + return self.__builtin_helper__( + self.get_local_data().__ifloordiv__, other, + inplace = True) def __mul__(self, other): return self.__builtin_helper__(self.get_local_data().__mul__, other) @@ -336,8 +360,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__rmul__, other) def __imul__(self, other): - return self.__inplace_builtin_helper__(self.get_local_data().__imul__, - other) + return self.__builtin_helper__(self.get_local_data().__imul__, + other, + inplace = True) def __pow__(self, other): return self.__builtin_helper__(self.get_local_data().__pow__, other) @@ -346,8 +371,9 @@ class distributed_data_object(object): return self.__builtin_helper__(self.get_local_data().__rpow__, other) def __ipow__(self, other): - return self.__inplace_builtin_helper__(self.get_local_data().__ipow__, - other) + return self.___builtin_helper__(self.get_local_data().__ipow__, + other, + inplace = True) def __len__(self): return self.shape[0] @@ -392,24 +418,30 @@ class distributed_data_object(object): def __setitem__(self, key, data): self.set_data(data, key) - def _minmaxhelper(self, function, **kwargs): + def _contraction_helper(self, function, **kwargs): local = function(self.data, **kwargs) local_list = self.distributor._allgather(local) global_ = function(local_list, axis=0) return global_ def amin(self, **kwargs): - return self._minmaxhelper(np.amin, **kwargs) + return self._contraction_helper(np.amin, **kwargs) def nanmin(self, **kwargs): - return self._minmaxhelper(np.nanmin, **kwargs) + return self._contraction_helper(np.nanmin, **kwargs) def amax(self, **kwargs): - return self._minmaxhelper(np.amax, **kwargs) + return self._contraction_helper(np.amax, **kwargs) def nanmax(self, **kwargs): - return self._minmaxhelper(np.nanmax, **kwargs) + return self._contraction_helper(np.nanmax, **kwargs) + def sum(self, **kwargs): + return self._contraction_helper(np.sum, **kwargs) + + def prod(self, **kwargs): + return self._contraction_helper(np.prod, **kwargs) + def mean(self, power=1): ## compute the local means and the weights for the mean-mean. local_mean = np.mean(self.data**power) @@ -731,8 +763,13 @@ class distributed_data_object(object): for i in sliceified: if i == True: temp_shape += (1,) + if data.shape[j] == 1: + j +=1 else: - temp_shape += (data.shape[j],) + try: + temp_shape += (data.shape[j],) + except(IndexError): + temp_shape += (1,) j += 1 ## take into account that the sliceified tuple may be too short, because ## of a non-exaustive list of slices diff --git a/nifty_tools.py b/nifty_tools.py deleted file mode 100644 index e66bee0d7fad1ae0a9990998941a29bfccb564fa..0000000000000000000000000000000000000000 --- a/nifty_tools.py +++ /dev/null @@ -1,650 +0,0 @@ -## 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/>. - -""" - .. __ ____ __ - .. /__/ / _/ / /_ - .. __ ___ __ / /_ / _/ __ __ - .. / _ | / / / _/ / / / / / / - .. / / / / / / / / / /_ / /_/ / - .. /__/ /__/ /__/ /__/ \___/ \___ / tools - .. /______/ - - This module extends NIFTY with a nifty set of tools including further - operators, namely the :py:class:`invertible_operator` and the - :py:class:`propagator_operator`, and minimization schemes, namely - :py:class:`steepest_descent` and :py:class:`conjugate_gradient`. Those - tools are supposed to support the user in solving information field - theoretical problems (almost) without numerical pain. - -""" -from __future__ import division -#from nifty_core import * -import numpy as np -from nifty_about import notification, about -from nifty_core import field -#from nifty_core import space, \ -# field -#from operators import operator, \ -# diagonal_operator - - - -##============================================================================= - -class conjugate_gradient(object): - """ - .. _______ ____ __ - .. / _____/ / _ / - .. / /____ __ / /_/ / __ - .. \______//__/ \____ //__/ class - .. /______/ - - NIFTY tool class for conjugate gradient - - This tool minimizes :math:`A x = b` with respect to `x` given `A` and - `b` using a conjugate gradient; i.e., a step-by-step minimization - relying on conjugated gradient directions. Further, `A` is assumed to - be a positive definite and self-adjoint operator. The use of a - preconditioner `W` that is roughly the inverse of `A` is optional. - For details on the methodology refer to [#]_, for details on usage and - output, see the notes below. - - Parameters - ---------- - A : {operator, function} - Operator `A` applicable to a field. - b : field - Resulting field of the operation `A(x)`. - W : {operator, function}, *optional* - Operator `W` that is a preconditioner on `A` and is applicable to a - field (default: None). - spam : function, *optional* - Callback function which is given the current `x` and iteration - counter each iteration (default: None). - reset : integer, *optional* - Number of iterations after which to restart; i.e., forget previous - conjugated directions (default: sqrt(b.dim())). - note : bool, *optional* - Indicates whether notes are printed or not (default: False). - - See Also - -------- - scipy.sparse.linalg.cg - - Notes - ----- - After initialization by `__init__`, the minimizer is started by calling - it using `__call__`, which takes additional parameters. Notifications, - if enabled, will state the iteration number, current step widths - `alpha` and `beta`, the current relative residual `delta` that is - compared to the tolerance, and the convergence level if changed. - The minimizer will exit in three states: DEAD if alpha becomes - infinite, QUIT if the maximum number of iterations is reached, or DONE - if convergence is achieved. Returned will be the latest `x` and the - latest convergence level, which can evaluate ``True`` for the exit - states QUIT and DONE. - - References - ---------- - .. [#] J. R. Shewchuk, 1994, `"An Introduction to the Conjugate - Gradient Method Without the Agonizing Pain" - <http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf>`_ - - Examples - -------- - >>> b = field(point_space(2), val=[1, 9]) - >>> A = diagonal_operator(b.domain, diag=[4, 3]) - >>> x,convergence = conjugate_gradient(A, b, note=True)(tol=1E-4, clevel=3) - iteration : 00000001 alpha = 3.3E-01 beta = 1.3E-03 delta = 3.6E-02 - iteration : 00000002 alpha = 2.5E-01 beta = 7.6E-04 delta = 1.0E-03 - iteration : 00000003 alpha = 3.3E-01 beta = 2.5E-04 delta = 1.6E-05 convergence level : 1 - iteration : 00000004 alpha = 2.5E-01 beta = 1.8E-06 delta = 2.1E-08 convergence level : 2 - iteration : 00000005 alpha = 2.5E-01 beta = 2.2E-03 delta = 1.0E-09 convergence level : 3 - ... done. - >>> bool(convergence) - True - >>> x.val # yields 1/4 and 9/3 - array([ 0.25, 3. ]) - - Attributes - ---------- - A : {operator, function} - Operator `A` applicable to a field. - x : field - Current field. - b : field - Resulting field of the operation `A(x)`. - W : {operator, function} - Operator `W` that is a preconditioner on `A` and is applicable to a - field; can be ``None``. - spam : function - Callback function which is given the current `x` and iteration - counter each iteration; can be ``None``. - reset : integer - Number of iterations after which to restart; i.e., forget previous - conjugated directions (default: sqrt(b.dim())). - note : notification - Notification instance. - - """ - def __init__(self,A,b,W=None,spam=None,reset=None,note=False): - """ - Initializes the conjugate_gradient and sets the attributes (except - for `x`). - - Parameters - ---------- - A : {operator, function} - Operator `A` applicable to a field. - b : field - Resulting field of the operation `A(x)`. - W : {operator, function}, *optional* - Operator `W` that is a preconditioner on `A` and is applicable to a - field (default: None). - spam : function, *optional* - Callback function which is given the current `x` and iteration - counter each iteration (default: None). - reset : integer, *optional* - Number of iterations after which to restart; i.e., forget previous - conjugated directions (default: sqrt(b.dim())). - note : bool, *optional* - Indicates whether notes are printed or not (default: False). - - """ - if(hasattr(A,"__call__")): - self.A = A ## applies A - else: - raise AttributeError(about._errors.cstring("ERROR: invalid input.")) - self.b = b - - if(W is None)or(hasattr(W,"__call__")): - self.W = W ## applies W ~ A_inverse - else: - raise AttributeError(about._errors.cstring("ERROR: invalid input.")) - - self.spam = spam ## serves as callback given x and iteration number - if(reset is None): ## 2 < reset ~ sqrt(dim) - self.reset = max(2,int(np.sqrt(b.domain.dim(split=False)))) - else: - self.reset = max(2,int(reset)) - self.note = notification(default=bool(note)) - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def __call__(self,x0=None,**kwargs): ## > runs cg with/without preconditioner - """ - Runs the conjugate gradient minimization. - - Parameters - ---------- - x0 : field, *optional* - Starting guess for the minimization. - tol : scalar, *optional* - Tolerance specifying convergence; measured by current relative - residual (default: 1E-4). - clevel : integer, *optional* - Number of times the tolerance should be undershot before - exiting (default: 1). - limii : integer, *optional* - Maximum number of iterations performed (default: 10 * b.dim()). - - Returns - ------- - x : field - Latest `x` of the minimization. - convergence : integer - Latest convergence level indicating whether the minimization - has converged or not. - - """ - self.x = field(self.b.domain,val=x0,target=self.b.target) - - if(self.W is None): - return self._calc_without(**kwargs) - else: - return self._calc_with(**kwargs) - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def _calc_without(self,tol=1E-4,clevel=1,limii=None): ## > runs cg without preconditioner - clevel = int(clevel) - if(limii is None): - limii = 10*self.b.domain.dim(split=False) - else: - limii = int(limii) - - r = self.b-self.A(self.x) - d = field(self.b.domain,val=np.copy(r.val),target=self.b.target) - gamma = r.dot(d) - if(gamma==0): - return self.x,clevel+1 - delta_ = np.absolute(gamma)**(-0.5) - - convergence = 0 - ii = 1 - while(True): - q = self.A(d) - alpha = gamma/d.dot(q) ## positive definite - if(not np.isfinite(alpha)): - self.note.cprint("\niteration : %08u alpha = NAN\n... dead."%ii) - return self.x,0 - self.x += alpha*d - if(np.signbit(np.real(alpha))): - about.warnings.cprint("WARNING: positive definiteness of A violated.") - r = self.b-self.A(self.x) - elif(ii%self.reset==0): - r = self.b-self.A(self.x) - else: - r -= alpha*q - gamma_ = gamma - gamma = r.dot(r) - beta = max(0,gamma/gamma_) ## positive definite - d = r+beta*d - - delta = delta_*np.absolute(gamma)**0.5 - self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta))) - if(gamma==0): - convergence = clevel+1 - self.note.cprint(" convergence level : INF\n... done.") - break - elif(np.absolute(delta)<tol): - convergence += 1 - self.note.cflush(" convergence level : %u"%convergence) - if(convergence==clevel): - self.note.cprint("\n... done.") - break - else: - convergence = max(0,convergence-1) - if(ii==limii): - self.note.cprint("\n... quit.") - break - - if(self.spam is not None): - self.spam(self.x,ii) - - ii += 1 - - if(self.spam is not None): - self.spam(self.x,ii) - - return self.x,convergence - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def _calc_with(self,tol=1E-4,clevel=1,limii=None): ## > runs cg with preconditioner - - clevel = int(clevel) - if(limii is None): - limii = 10*self.b.domain.dim(split=False) - else: - limii = int(limii) - - r = self.b-self.A(self.x) - d = self.W(r) - gamma = r.dot(d) - if(gamma==0): - return self.x,clevel+1 - delta_ = np.absolute(gamma)**(-0.5) - - convergence = 0 - ii = 1 - while(True): - q = self.A(d) - alpha = gamma/d.dot(q) ## positive definite - if(not np.isfinite(alpha)): - self.note.cprint("\niteration : %08u alpha = NAN\n... dead."%ii) - return self.x,0 - self.x += alpha*d ## update - if(np.signbit(np.real(alpha))): - about.warnings.cprint("WARNING: positive definiteness of A violated.") - r = self.b-self.A(self.x) - elif(ii%self.reset==0): - r = self.b-self.A(self.x) - else: - r -= alpha*q - s = self.W(r) - gamma_ = gamma - gamma = r.dot(s) - if(np.signbit(np.real(gamma))): - about.warnings.cprint("WARNING: positive definiteness of W violated.") - beta = max(0,gamma/gamma_) ## positive definite - d = s+beta*d ## conjugated gradient - - delta = delta_*np.absolute(gamma)**0.5 - self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta))) - if(gamma==0): - convergence = clevel+1 - self.note.cprint(" convergence level : INF\n... done.") - break - elif(np.absolute(delta)<tol): - convergence += 1 - self.note.cflush(" convergence level : %u"%convergence) - if(convergence==clevel): - self.note.cprint("\n... done.") - break - else: - convergence = max(0,convergence-1) - if(ii==limii): - self.note.cprint("\n... quit.") - break - - if(self.spam is not None): - self.spam(self.x,ii) - - ii += 1 - - if(self.spam is not None): - self.spam(self.x,ii) - - return self.x,convergence - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def __repr__(self): - return "<nifty_tools.conjugate_gradient>" - -##============================================================================= - - - - - -##============================================================================= - -class steepest_descent(object): - """ - .. __ - .. / / - .. _______ ____/ / - .. / _____/ / _ / - .. /_____ / __ / /_/ / __ - .. /_______//__/ \______|/__/ class - - NIFTY tool class for steepest descent minimization - - This tool minimizes a scalar energy-function by steepest descent using - the functions gradient. Steps and step widths are choosen according to - the Wolfe conditions [#]_. For details on usage and output, see the - notes below. - - Parameters - ---------- - eggs : function - Given the current `x` it returns the tuple of energy and gradient. - spam : function, *optional* - Callback function which is given the current `x` and iteration - counter each iteration (default: None). - a : {4-tuple}, *optional* - Numbers obeying 0 < a1 ~ a2 < 1 ~ a3 < a4 that modify the step - widths (default: (0.2,0.5,1,2)). - c : {2-tuple}, *optional* - Numbers obeying 0 < c1 < c2 < 1 that specify the Wolfe-conditions - (default: (1E-4,0.9)). - note : bool, *optional* - Indicates whether notes are printed or not (default: False). - - See Also - -------- - scipy.optimize.fmin_cg, scipy.optimize.fmin_ncg, - scipy.optimize.fmin_l_bfgs_b - - Notes - ----- - After initialization by `__init__`, the minimizer is started by calling - it using `__call__`, which takes additional parameters. Notifications, - if enabled, will state the iteration number, current step width `alpha`, - current maximal change `delta` that is compared to the tolerance, and - the convergence level if changed. The minimizer will exit in three - states: DEAD if no step width above 1E-13 is accepted, QUIT if the - maximum number of iterations is reached, or DONE if convergence is - achieved. Returned will be the latest `x` and the latest convergence - level, which can evaluate ``True`` for all exit states. - - References - ---------- - .. [#] J. Nocedal and S. J. Wright, Springer 2006, "Numerical - Optimization", ISBN: 978-0-387-30303-1 (print) / 978-0-387-40065-5 - `(online) <http://link.springer.com/book/10.1007/978-0-387-40065-5/page/1>`_ - - Examples - -------- - >>> def egg(x): - ... E = 0.5*x.dot(x) # energy E(x) -- a two-dimensional parabola - ... g = x # gradient - ... return E,g - >>> x = field(point_space(2), val=[1, 3]) - >>> x,convergence = steepest_descent(egg, note=True)(x0=x, tol=1E-4, clevel=3) - iteration : 00000001 alpha = 1.0E+00 delta = 6.5E-01 - iteration : 00000002 alpha = 2.0E+00 delta = 1.4E-01 - iteration : 00000003 alpha = 1.6E-01 delta = 2.1E-03 - iteration : 00000004 alpha = 2.6E-03 delta = 3.0E-04 - iteration : 00000005 alpha = 2.0E-04 delta = 5.3E-05 convergence level : 1 - iteration : 00000006 alpha = 8.2E-05 delta = 4.4E-06 convergence level : 2 - iteration : 00000007 alpha = 6.6E-06 delta = 3.1E-06 convergence level : 3 - ... done. - >>> bool(convergence) - True - >>> x.val # approximately zero - array([ -6.87299426e-07 -2.06189828e-06]) - - Attributes - ---------- - x : field - Current field. - eggs : function - Given the current `x` it returns the tuple of energy and gradient. - spam : function - Callback function which is given the current `x` and iteration - counter each iteration; can be ``None``. - a : {4-tuple} - Numbers obeying 0 < a1 ~ a2 < 1 ~ a3 < a4 that modify the step - widths (default: (0.2,0.5,1,2)). - c : {2-tuple} - Numbers obeying 0 < c1 < c2 < 1 that specify the Wolfe-conditions - (default: (1E-4,0.9)). - note : notification - Notification instance. - - """ - def __init__(self,eggs,spam=None,a=(0.2,0.5,1,2),c=(1E-4,0.9),note=False): - """ - Initializes the steepest_descent and sets the attributes (except - for `x`). - - Parameters - ---------- - eggs : function - Given the current `x` it returns the tuple of energy and gradient. - spam : function, *optional* - Callback function which is given the current `x` and iteration - counter each iteration (default: None). - a : {4-tuple}, *optional* - Numbers obeying 0 < a1 ~ a2 < 1 ~ a3 < a4 that modify the step - widths (default: (0.2,0.5,1,2)). - c : {2-tuple}, *optional* - Numbers obeying 0 < c1 < c2 < 1 that specify the Wolfe-conditions - (default: (1E-4,0.9)). - note : bool, *optional* - Indicates whether notes are printed or not (default: False). - - """ - self.eggs = eggs ## returns energy and gradient - - self.spam = spam ## serves as callback given x and iteration number - self.a = a ## 0 < a1 ~ a2 < 1 ~ a3 < a4 - self.c = c ## 0 < c1 < c2 < 1 - self.note = notification(default=bool(note)) - - self._alpha = None ## last alpha - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def __call__(self,x0,alpha=1,tol=1E-4,clevel=8,limii=100000): - """ - Runs the steepest descent minimization. - - Parameters - ---------- - x0 : field - Starting guess for the minimization. - alpha : scalar, *optional* - Starting step width to be multiplied with normalized gradient - (default: 1). - tol : scalar, *optional* - Tolerance specifying convergence; measured by maximal change in - `x` (default: 1E-4). - clevel : integer, *optional* - Number of times the tolerance should be undershot before - exiting (default: 8). - limii : integer, *optional* - Maximum number of iterations performed (default: 100,000). - - Returns - ------- - x : field - Latest `x` of the minimization. - convergence : integer - Latest convergence level indicating whether the minimization - has converged or not. - - """ - if(not isinstance(x0,field)): - raise TypeError(about._errors.cstring("ERROR: invalid input.")) - self.x = x0 - - ## check for exsisting alpha - if(alpha is None): - if(self._alpha is not None): - alpha = self._alpha - else: - alpha = 1 - - clevel = max(1,int(clevel)) - limii = int(limii) - - E,g = self.eggs(self.x) ## energy and gradient - norm = g.norm() ## gradient norm - if(norm==0): - self.note.cprint("\niteration : 00000000 alpha = 0.0E+00 delta = 0.0E+00\n... done.") - return self.x,clevel+2 - - convergence = 0 - ii = 1 - while(True): - x_,E,g,alpha,a = self._get_alpha(E,g,norm,alpha) ## "news",alpha,a - - if(alpha is None): - self.note.cprint("\niteration : %08u alpha < 1.0E-13\n... dead."%ii) - break - else: - delta = np.absolute(g.val).max()*(alpha/norm) - self.note.cflush("\niteration : %08u alpha = %3.1E delta = %3.1E"%(ii,alpha,delta)) - ## update - self.x = x_ - alpha *= a - - norm = g.norm() ## gradient norm - if(delta==0): - convergence = clevel+2 - self.note.cprint(" convergence level : %u\n... done."%convergence) - break - elif(delta<tol): - convergence += 1 - self.note.cflush(" convergence level : %u"%convergence) - if(convergence==clevel): - convergence += int(ii==clevel) - self.note.cprint("\n... done.") - break - else: - convergence = max(0,convergence-1) - if(ii==limii): - self.note.cprint("\n... quit.") - break - - if(self.spam is not None): - self.spam(self.x,ii) - - ii += 1 - - if(self.spam is not None): - self.spam(self.x,ii) - - ## memorise last alpha - if(alpha is not None): - self._alpha = alpha/a ## undo update - - return self.x,convergence - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def _get_alpha(self,E,g,norm,alpha): ## > determines the new alpha - while(True): - ## Wolfe conditions - wolfe,x_,E_,g_,a = self._check_wolfe(E,g,norm,alpha) -# wolfe,x_,E_,g_,a = self._check_strong_wolfe(E,g,norm,alpha) - if(wolfe): - return x_,E_,g_,alpha,a - else: - alpha *= a - if(alpha<1E-13): - return None,None,None,None,None - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def _check_wolfe(self,E,g,norm,alpha): ## > checks the Wolfe conditions - x_ = self._get_x(g,norm,alpha) - pg = norm - E_,g_ = self.eggs(x_) - if(E_>E+self.c[0]*alpha*pg): - if(E_<E): - return True,x_,E_,g_,self.a[1] - return False,None,None,None,self.a[0] - pg_ = g.dot(g_)/norm - if(pg_<self.c[1]*pg): - return True,x_,E_,g_,self.a[3] - return True,x_,E_,g_,self.a[2] - -# def _check_strong_wolfe(self,E,g,norm,alpha): ## > checks the strong Wolfe conditions -# x_ = self._get_x(g,norm,alpha) -# pg = norm -# E_,g_ = self.eggs(x_) -# if(E_>E+self.c[0]*alpha*pg): -# if(E_<E): -# return True,x_,E_,g_,self.a[1] -# return False,None,None,None,self.a[0] -# apg_ = np.absolute(g.dot(g_))/norm -# if(apg_>self.c[1]*np.absolute(pg)): -# return True,x_,E_,g_,self.a[3] -# return True,x_,E_,g_,self.a[2] - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def _get_x(self,g,norm,alpha): ## > updates x - return self.x-g*(alpha/norm) - - ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - def __repr__(self): - return "<nifty_tools.steepest_descent>" - -##============================================================================= - diff --git a/operators/nifty_operators.py b/operators/nifty_operators.py index 894466601259186e243a4278725a8f6cee4d64cf..35f1269baee4e3fb9046db3f95aaa25b7d1bfb2c 100644 --- a/operators/nifty_operators.py +++ b/operators/nifty_operators.py @@ -26,7 +26,7 @@ from nifty.nifty_core import space, \ point_space, \ nested_space, \ field -from nifty.nifty_tools import conjugate_gradient +from nifty_minimization import conjugate_gradient from nifty_probing import trace_probing, \ diagonal_probing @@ -1803,8 +1803,9 @@ class power_operator(diagonal_operator): try: #diag = self.domain.enforce_power(spec,size=np.max(pindex,axis=None,out=None)+1)[pindex] temp_spec = self.domain.enforce_power( - spec,size=np.max(pindex,axis=None,out=None)+1) - diag = pindex.apply_scalar_function(lambda x: temp_spec[x]) + spec,size=np.max(pindex,axis=None,out=None)+1) + diag = pindex.apply_scalar_function(lambda x: temp_spec[x], + dtype = temp_spec.dtype.type) except(AttributeError): raise ValueError(about._errors.cstring("ERROR: invalid input.")) ## weight if ... diff --git a/rg/fft_rg.py b/rg/fft_rg.py index c092c55afe8d676db14ab4bdfe90dab6a15bd89b..cf626b963d6dcef987cbf92f273e4e52393151c2 100644 --- a/rg/fft_rg.py +++ b/rg/fft_rg.py @@ -8,7 +8,7 @@ from nifty.nifty_about import about # If this fails fall back to local gfft_rg try: - import pyfftw_BAD + import pyfftw fft_machine='pyfftw' except(ImportError): try: @@ -182,6 +182,10 @@ if fft_machine == 'pyfftw': ## cast input to_center = np.array(to_center_input) dimensions = np.array(dimensions_input) + + ## if none of the dimensions are zero centered, return a 1 + if np.all(to_center == 0): + return 1 if np.all(dimensions == np.array(1)) or \ np.all(dimensions == np.array([1])): @@ -221,7 +225,9 @@ if fft_machine == 'pyfftw': offset.reshape(offset.shape + \ (1,)*(np.array(args).ndim - 1)),1)),\ (2,)*to_center.size) - + ## Cast the core to the smallest integers we can get + core = core.astype(np.int8) + centering_mask = np.tile(core,dimensions//2) ## for the dimensions of odd size corresponding slices must be added for i in range(centering_mask.ndim): @@ -249,7 +255,7 @@ if fft_machine == 'pyfftw': def _get_plan_and_info(self,domain,codomain,**kwargs): ## generate a id-tuple which identifies the domain-codomain setting - temp_id = (domain.__identifier__(), codomain.__identifier__()) + temp_id = (domain._identifier(), codomain._identifier()) ## generate the plan_and_info object if not already there if not temp_id in self.plan_dict: self.plan_dict[temp_id]=_fftw_plan_and_info(domain, codomain, diff --git a/rg/nifty_power_conversion_rg.py b/rg/nifty_power_conversion_rg.py index d478ffca0cf88da86212ab10f5f6fe3f635c6fa7..9d8c6c60ea8a95e2f1932f6695f41800140f4ef2 100644 --- a/rg/nifty_power_conversion_rg.py +++ b/rg/nifty_power_conversion_rg.py @@ -106,7 +106,7 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True): return logmean.real,p1.real -def power_forward_conversion_rg(k_space,p,mean=0,bare=True): +def power_forward_conversion_rg(k_space, p, mean=0, bare=True): """ This function is designed to convert a theoretical/statistical power spectrum of a Gaussian field to the theoretical power spectrum of @@ -137,6 +137,44 @@ def power_forward_conversion_rg(k_space,p,mean=0,bare=True): `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_ """ + pindex = k_space.power_indices['pindex'] + weight = k_space.get_weight() + ## Cast the supplied spectrum + spec = k_space.enforce_power(p) + ## Now we mimick the weightning behaviour of + ## spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False) + ## by appliying the weight from the k_space + if bare == True: + spec *= weight + + S_val = pindex.apply_scalar_function(lambda x: spec[x], + dtype=spec.dtype.type) + + ## S_x is a field + S_x = field(k_space, val=S_val, zerocenter=True).transform() + ## s_0 is a scalar + s_0 = k_space.calc_weight(S_val, power = 1).sum() + + ## Calculate the new power_field + S_x += s_0 + S_x += 2*mean + + print S_x + print s_0 + power_field = S_x.apply_scalar_function(np.exp, inplace=True) + + new_spec = power_field.power()**(0.5) + + ## Mimik + ## power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real + if bare == True: + new_spec /= weight + + return new_spec.real + + + """ + pindex = k_space.get_power_indices()[2] spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False) @@ -154,3 +192,4 @@ def power_forward_conversion_rg(k_space,p,mean=0,bare=True): else: return p1.real + """ \ No newline at end of file diff --git a/rg/nifty_rg.py b/rg/nifty_rg.py index f1551441a8ec7de7ca019b4fb33002e3a01327f4..499c285551aa1b7b4839d732b96d09543989663d 100644 --- a/rg/nifty_rg.py +++ b/rg/nifty_rg.py @@ -583,9 +583,9 @@ class rg_space(point_space): "ERROR: Data has incompatible shape!")) ## Check the datatype - if x.dtype != self.datatype: + if x.dtype < self.datatype: about.warnings.cflush(\ - "WARNING: Datatypes are uneqal (own: "\ + "WARNING: Datatypes are uneqal/of conflicting precision (own: "\ + str(self.datatype) + " <> foreign: " + str(x.dtype) \ + ") and will be casted! "\ + "Potential loss of precision!\n") @@ -1161,13 +1161,18 @@ class rg_space(point_space): y : numpy.ndarray Weighted array. """ - x = self.cast(x) - is_hermitianQ = x.hermitian + #x = self.cast(x) + if isinstance(x, distributed_data_object): + is_hermitianQ = x.hermitian ## weight - x = x * (np.prod(self.vol)**power) - x.hermitian = is_hermitianQ + x = x * self.get_weight(power = power) + if isinstance(x, distributed_data_object): + x.hermitian = is_hermitianQ return x + def get_weight(self, power = 1): + return np.prod(self.vol)**power + ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def calc_dot(self, x, y): """ @@ -1426,11 +1431,7 @@ class rg_space(point_space): fieldabs = abs(x)**2 power_spectrum = np.zeros(rho.shape) - """ - ##TODO: Replace this super slow ndindex solution - for ii in np.ndindex(pindex.shape): - power_spectrum[pindex[ii]] += fieldabs[ii] - """ + ## In order to make the summation over identical pindices fast, ## the pindex and the kindex must have the same distribution strategy @@ -1446,20 +1447,7 @@ class rg_space(point_space): power_spectrum =\ pindex.distributor._allgather(local_power_spectrum) power_spectrum = np.sum(power_spectrum, axis = 0) - - """ - ## Iterate over the k-vectors, extract those fieldabs, where the pindex - ## has the according value and build the sum of the resulting array - - power_spectrum = np.zeros(rho.size, dtype = np.float) - - - for ii in xrange(rho.size): - ## extract those fieldabs where the pindex equals the current ii - extracted_fieldabs = working_field[pindex == ii] - ## sum the extracted field values up and store them - power_spectrum[ii] = np.sum(extracted_fieldabs) - """ + ## Divide out the degeneracy factor power_spectrum /= rho return power_spectrum diff --git a/rg/powerspectrum.py b/rg/powerspectrum.py deleted file mode 100644 index 6d2c812be3abe69dfd33bcace92788b626ba8493..0000000000000000000000000000000000000000 --- a/rg/powerspectrum.py +++ /dev/null @@ -1,761 +0,0 @@ -## 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/>. - - -from __future__ import division -import numpy as np - - - - -def draw_vector_nd(axes,dgrid,ps,symtype=0,fourier=False,zerocentered=False,kpack=None): - - """ - Draws a n-dimensional field on a regular grid from a given power - spectrum. The grid parameters need to be specified, together with a - couple of global options explained below. The dimensionality of the - field is determined automatically. - - Parameters - ---------- - axes : ndarray - An array with the length of each axis. - - dgrid : ndarray - An array with the pixel length of each axis. - - ps : ndarray - The power spectrum as a function of Fourier modes. - - symtype : int {0,1,2} : *optional* - Whether the output should be real valued (0), complex-hermitian (1) - or complex without symmetry (2). (default=0) - - fourier : bool : *optional* - Whether the output should be in Fourier space or not - (default=False). - - zerocentered : bool : *optional* - Whether the output array should be zerocentered, i.e. starting with - negative Fourier modes going over the zero mode to positive modes, - or not zerocentered, where zero, positive and negative modes are - simpy ordered consecutively. - - Returns - ------- - field : ndarray - The drawn random field. - - """ - if(kpack is None): - kdict = np.fft.fftshift(nkdict_fast(axes,dgrid,fourier)) - klength = nklength(kdict) - else: - kdict = kpack[1][np.fft.ifftshift(kpack[0],axes=shiftaxes(zerocentered,st_to_zero_mode=False))] - klength = kpack[1] - - #output is in position space - if(not fourier): - - #output is real-valued - if(symtype==0): - vector = drawherm(klength,kdict,ps) - if(np.any(zerocentered==True)): - return np.real(np.fft.fftshift(np.fft.ifftn(vector),axes=shiftaxes(zerocentered))) - else: - return np.real(np.fft.ifftn(vector)) - - #output is complex with hermitian symmetry - elif(symtype==1): - vector = drawwild(klength,kdict,ps,real_corr=2) - if(np.any(zerocentered==True)): - return np.fft.fftshift(np.fft.ifftn(np.real(vector)),axes=shiftaxes(zerocentered)) - else: - return np.fft.ifftn(np.real(vector)) - - #output is complex without symmetry - else: - vector = drawwild(klength,kdict,ps) - if(np.any(zerocentered==True)): - return np.fft.fftshift(np.fft.ifftn(vector),axes=shiftaxes(zerocentered)) - else: - return np.fft.ifftn(vector) - - #output is in fourier space - else: - - #output is real-valued - if(symtype==0): - vector = drawwild(klength,kdict,ps,real_corr=2) - if np.any(zerocentered == True): - return np.real(np.fft.fftshift(vector,axes=shiftaxes(zerocentered))) - else: - return np.real(vector) - - #output is complex with hermitian symmetry - elif(symtype==1): - vector = drawherm(klength,kdict,ps) - if(np.any(zerocentered==True)): - return np.fft.fftshift(vector,axes=shiftaxes(zerocentered)) - else: - return vector - - #output is complex without symmetry - else: - vector = drawwild(klength,kdict,ps) - if(np.any(zerocentered==True)): - return np.fft.fftshift(vector,axes=shiftaxes(zerocentered)) - else: - return vector - - -#def calc_ps(field,axes,dgrid,zerocentered=False,fourier=False): -# -# """ -# Calculates the power spectrum of a given field assuming that the field -# is statistically homogenous and isotropic. -# -# Parameters -# ---------- -# field : ndarray -# The input field from which the power spectrum should be determined. -# -# axes : ndarray -# An array with the length of each axis. -# -# dgrid : ndarray -# An array with the pixel length of each axis. -# -# zerocentered : bool : *optional* -# Whether the output array should be zerocentered, i.e. starting with -# negative Fourier modes going over the zero mode to positive modes, -# or not zerocentered, where zero, positive and negative modes are -# simpy ordered consecutively. -# -# fourier : bool : *optional* -# Whether the output should be in Fourier space or not -# (default=False). -# -# """ -# -# ## field absolutes -# if(not fourier): -# foufield = np.fft.fftshift(np.fft.fftn(field)) -# elif(np.any(zerocentered==False)): -# foufield = np.fft.fftshift(field, axes=shiftaxes(zerocentered,st_to_zero_mode=True)) -# else: -# foufield = field -# fieldabs = np.abs(foufield)**2 -# -# kdict = nkdict_fast(axes,dgrid,fourier) -# klength = nklength(kdict) -# -# ## power spectrum -# ps = np.zeros(klength.size) -# rho = np.zeros(klength.size) -# for ii in np.ndindex(kdict.shape): -# position = np.searchsorted(klength,kdict[ii]) -# rho[position] += 1 -# ps[position] += fieldabs[ii] -# ps = np.divide(ps,rho) -# return ps - -def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,kindex=None,rho=None): - - """ - Calculates the power spectrum of a given field faster assuming that the - field is statistically homogenous and isotropic. - - Parameters - ---------- - field : ndarray - The input field from which the power spectrum should be determined. - - axes : ndarray - An array with the length of each axis. - - dgrid : ndarray - An array with the pixel length of each axis. - - zerocentered : bool : *optional* - Whether the output array should be zerocentered, i.e. starting with - negative Fourier modes going over the zero mode to positive modes, - or not zerocentered, where zero, positive and negative modes are - simpy ordered consecutively. - - fourier : bool : *optional* - Whether the output should be in Fourier space or not - (default=False). - - pindex : ndarray - Index of the Fourier grid points in a numpy.ndarray ordered - following the zerocentered flag (default=None). - - kindex : ndarray - Array of all k-vector lengths (default=None). - - rho : ndarray - Degeneracy of the Fourier grid, indicating how many k-vectors in - Fourier space have the same length (default=None). - - """ - ## field absolutes - if(not fourier): - foufield = np.fft.fftshift(np.fft.fftn(field)) - elif(np.any(zerocentered==False)): - foufield = np.fft.fftshift(field, axes=shiftaxes(zerocentered,st_to_zero_mode=True)) - else: - foufield = field - fieldabs = np.abs(foufield)**2 - - if(rho is None): - if(pindex is None): - ## kdict - kdict = nkdict_fast(axes,dgrid,fourier) - ## klength - if(kindex is None): - klength = nklength(kdict) - else: - klength = kindex - ## power spectrum - ps = np.zeros(klength.size) - rho = np.zeros(klength.size) - for ii in np.ndindex(kdict.shape): - position = np.searchsorted(klength,kdict[ii]) - ps[position] += fieldabs[ii] - rho[position] += 1 - else: - ## zerocenter pindex - if(np.any(zerocentered==False)): - pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True)) - ## power spectrum - ps = np.zeros(np.max(pindex)+1) - rho = np.zeros(ps.size) - for ii in np.ndindex(pindex.shape): - ps[pindex[ii]] += fieldabs[ii] - rho[pindex[ii]] += 1 - elif(pindex is None): - ## kdict - kdict = nkdict_fast(axes,dgrid,fourier) - ## klength - if(kindex is None): - klength = nklength(kdict) - else: - klength = kindex - ## power spectrum - ps = np.zeros(klength.size) - for ii in np.ndindex(kdict.shape): - position = np.searchsorted(klength,kdict[ii]) - ps[position] += fieldabs[ii] - else: - ## zerocenter pindex - if(np.any(zerocentered==False)): - pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True)) - ## power spectrum - ps = np.zeros(rho.size) - for ii in np.ndindex(pindex.shape): - ps[pindex[ii]] += fieldabs[ii] - - ps = np.divide(ps,rho) - return ps - -# -#def get_power_index(axes,dgrid,zerocentered,irred=False,fourier=True): -# -# """ -# Returns the index of the Fourier grid points in a numpy -# array, ordered following the zerocentered flag. -# -# Parameters -# ---------- -# axes : ndarray -# An array with the length of each axis. -# -# dgrid : ndarray -# An array with the pixel length of each axis. -# -# zerocentered : bool -# Whether the output array should be zerocentered, i.e. starting with -# negative Fourier modes going over the zero mode to positive modes, -# or not zerocentered, where zero, positive and negative modes are -# simpy ordered consecutively. -# -# irred : bool : *optional* -# If True, the function returns an array of all k-vector lengths and -# their degeneracy factors. If False, just the power index array is -# returned. -# -# fourier : bool : *optional* -# Whether the output should be in Fourier space or not -# (default=False). -# -# Returns -# ------- -# index or {klength, rho} : scalar or list -# Returns either an array of all k-vector lengths and -# their degeneracy factors or just the power index array -# depending on the flag irred. -# -# """ -# -# ## kdict, klength -# if(np.any(zerocentered==False)): -# kdict = np.fft.fftshift(nkdict_fast(axes,dgrid,fourier),axes=shiftaxes(zerocentered,st_to_zero_mode=True)) -# else: -# kdict = nkdict_fast(axes,dgrid,fourier) -# klength = nklength(kdict) -# ## output -# if(irred): -# rho = np.zeros(klength.shape,dtype=np.int) -# for ii in np.ndindex(kdict.shape): -# rho[np.searchsorted(klength,kdict[ii])] += 1 -# return klength,rho -# else: -# ind = np.empty(axes,dtype=np.int) -# for ii in np.ndindex(kdict.shape): -# ind[ii] = np.searchsorted(klength,kdict[ii]) -# return ind -# -# -#def get_power_indices(axes,dgrid,zerocentered,fourier=True): -# """ -# Returns the index of the Fourier grid points in a numpy -# array, ordered following the zerocentered flag. -# -# Parameters -# ---------- -# axes : ndarray -# An array with the length of each axis. -# -# dgrid : ndarray -# An array with the pixel length of each axis. -# -# zerocentered : bool -# Whether the output array should be zerocentered, i.e. starting with -# negative Fourier modes going over the zero mode to positive modes, -# or not zerocentered, where zero, positive and negative modes are -# simpy ordered consecutively. -# -# irred : bool : *optional* -# If True, the function returns an array of all k-vector lengths and -# their degeneracy factors. If False, just the power index array is -# returned. -# -# fourier : bool : *optional* -# Whether the output should be in Fourier space or not -# (default=False). -# -# Returns -# ------- -# index, klength, rho : ndarrays -# Returns the power index array, an array of all k-vector lengths and -# their degeneracy factors. -# -# """ -# -# ## kdict, klength -# if(np.any(zerocentered==False)): -# kdict = np.fft.fftshift(nkdict_fast(axes,dgrid,fourier),axes=shiftaxes(zerocentered,st_to_zero_mode=True)) -# else: -# kdict = nkdict_fast(axes,dgrid,fourier) -# klength = nklength(kdict) -# ## output -# ind = np.empty(axes,dtype=np.int) -# rho = np.zeros(klength.shape,dtype=np.int) -# for ii in np.ndindex(kdict.shape): -# ind[ii] = np.searchsorted(klength,kdict[ii]) -# rho[ind[ii]] += 1 -# return ind,klength,rho -# -# -#def get_power_indices2(axes,dgrid,zerocentered,fourier=True): -# """ -# Returns the index of the Fourier grid points in a numpy -# array, ordered following the zerocentered flag. -# -# Parameters -# ---------- -# axes : ndarray -# An array with the length of each axis. -# -# dgrid : ndarray -# An array with the pixel length of each axis. -# -# zerocentered : bool -# Whether the output array should be zerocentered, i.e. starting with -# negative Fourier modes going over the zero mode to positive modes, -# or not zerocentered, where zero, positive and negative modes are -# simpy ordered consecutively. -# -# irred : bool : *optional* -# If True, the function returns an array of all k-vector lengths and -# their degeneracy factors. If False, just the power index array is -# returned. -# -# fourier : bool : *optional* -# Whether the output should be in Fourier space or not -# (default=False). -# -# Returns -# ------- -# index, klength, rho : ndarrays -# Returns the power index array, an array of all k-vector lengths and -# their degeneracy factors. -# -# """ -# -# ## kdict, klength -# if(np.any(zerocentered==False)): -# kdict = np.fft.fftshift(nkdict_fast2(axes,dgrid,fourier),axes=shiftaxes(zerocentered,st_to_zero_mode=True)) -# else: -# kdict = nkdict_fast2(axes,dgrid,fourier) -# -# klength,rho,ind = nkdict_to_indices(kdict) -# -# return ind,klength,rho -# -#def nkdict_to_indices(kdict): -# -# kindex,pindex = np.unique(kdict,return_inverse=True) -# pindex = pindex.reshape(kdict.shape) -# -# rho = pindex.flatten() -# rho.sort() -# rho = np.unique(rho,return_index=True,return_inverse=False)[1] -# rho = np.append(rho[1:]-rho[:-1],[np.prod(pindex.shape)-rho[-1]]) -# -# return kindex,rho,pindex -# -# -# -#def bin_power_indices(pindex,kindex,rho,log=False,nbin=None,binbounds=None): -# """ -# Returns the (re)binned power indices associated with the Fourier grid. -# -# Parameters -# ---------- -# pindex : ndarray -# Index of the Fourier grid points in a numpy.ndarray ordered -# following the zerocentered flag (default=None). -# kindex : ndarray -# Array of all k-vector lengths (default=None). -# rho : ndarray -# Degeneracy of the Fourier grid, indicating how many k-vectors in -# Fourier space have the same length (default=None). -# log : bool -# Flag specifying if the binning is performed on logarithmic scale -# (default: False). -# nbin : integer -# Number of used bins (default: None). -# binbounds : {list, array} -# Array-like inner boundaries of the used bins (default: None). -# -# Returns -# ------- -# pindex, kindex, rho : ndarrays -# The (re)binned power indices. -# -# """ -# ## boundaries -# if(binbounds is not None): -# binbounds = np.sort(binbounds) -# ## equal binning -# else: -# if(log is None): -# log = False -# if(log): -# k = np.r_[0,np.log(kindex[1:])] -# else: -# k = kindex -# dk = np.max(k[2:]-k[1:-1]) ## minimal dk -# if(nbin is None): -# nbin = int((k[-1]-0.5*(k[2]+k[1]))/dk-0.5) ## maximal nbin -# else: -# nbin = min(int(nbin),int((k[-1]-0.5*(k[2]+k[1]))/dk+2.5)) -# dk = (k[-1]-0.5*(k[2]+k[1]))/(nbin-2.5) -# binbounds = np.r_[0.5*(3*k[1]-k[2]),0.5*(k[1]+k[2])+dk*np.arange(nbin-2)] -# if(log): -# binbounds = np.exp(binbounds) -# ## reordering -# reorder = np.searchsorted(binbounds,kindex) -# rho_ = np.zeros(len(binbounds)+1,dtype=rho.dtype) -# kindex_ = np.empty(len(binbounds)+1,dtype=kindex.dtype) -# for ii in range(len(reorder)): -# if(rho_[reorder[ii]]==0): -# kindex_[reorder[ii]] = kindex[ii] -# rho_[reorder[ii]] += rho[ii] -# else: -# kindex_[reorder[ii]] = (kindex_[reorder[ii]]*rho_[reorder[ii]]+kindex[ii]*rho[ii])/(rho_[reorder[ii]]+rho[ii]) -# rho_[reorder[ii]] += rho[ii] -# -# return reorder[pindex],kindex_,rho_ -# - - -def nhermitianize(field,zerocentered): - - """ - Hermitianizes an arbitrary n-dimensional field. Becomes relatively slow - for large n. - - Parameters - ---------- - field : ndarray - The input field that should be hermitianized. - - zerocentered : bool - Whether the output array should be zerocentered, i.e. starting with - negative Fourier modes going over the zero mode to positive modes, - or not zerocentered, where zero, positive and negative modes are - simpy ordered consecutively. - - Returns - ------- - hermfield : ndarray - The hermitianized field. - - """ - ## shift zerocentered axes - if(np.any(zerocentered==True)): - field = np.fft.fftshift(field, axes=shiftaxes(zerocentered)) -# for index in np.ndenumerate(field): -# negind = tuple(-np.array(index[0])) -# field[negind] = np.conjugate(index[1]) -# if(field[negind]==field[index[0]]): -# field[index[0]] = np.abs(index[1])*(np.sign(index[1].real)+(np.sign(index[1].real)==0)*np.sign(index[1].imag)).astype(np.int) - subshape = np.array(field.shape,dtype=np.int) ## == axes - maxindex = subshape//2 - subshape[np.argmax(subshape)] = subshape[np.argmax(subshape)]//2+1 ## ~half larges axis - for ii in np.ndindex(tuple(subshape)): - negii = tuple(-np.array(ii)) - field[negii] = np.conjugate(field[ii]) - for ii in np.ndindex((2,)*maxindex.size): - index = tuple(ii*maxindex) - field[index] = np.abs(field[index])*(np.sign(field[index].real)+(np.sign(field[index].real)==0)*-np.sign(field[index].imag)).astype(np.int) ## minus since overwritten before - ## reshift zerocentered axes - if(np.any(zerocentered==True)): - field = np.fft.fftshift(field,axes=shiftaxes(zerocentered)) - return field - -def nhermitianize_fast(field,zerocentered,special=False): - - """ - Hermitianizes an arbitrary n-dimensional field faster. - Still becomes comparably slow for large n. - - Parameters - ---------- - field : ndarray - The input field that should be hermitianized. - - zerocentered : bool - Whether the output array should be zerocentered, i.e. starting with - negative Fourier modes going over the zero mode to positive modes, - or not zerocentered, where zero, positive and negative modes are - simpy ordered consecutively. - - special : bool, *optional* - Must be True for random fields drawn from Gaussian or pm1 - distributions. - - Returns - ------- - hermfield : ndarray - The hermitianized field. - - """ - ## shift zerocentered axes - if(np.any(zerocentered==True)): - field = np.fft.fftshift(field, axes=shiftaxes(zerocentered)) - dummy = np.conjugate(field) - ## mirror conjugate field - for ii in range(field.ndim): - dummy = np.swapaxes(dummy,0,ii) - dummy = np.flipud(dummy) - dummy = np.roll(dummy,1,axis=0) - dummy = np.swapaxes(dummy,0,ii) - if(special): ## special normalisation for certain random fields - field = np.sqrt(0.5)*(field+dummy) - maxindex = np.array(field.shape,dtype=np.int)//2 - print ('maxindex: ', maxindex) - for ii in np.ndindex((2,)*maxindex.size): - print ('ii: ', ii) - index = tuple(ii*maxindex) - print ('index: ', index) - field[index] *= np.sqrt(0.5) - else: ## regular case - field = 0.5*(field+dummy) - #field = dummy - ## reshift zerocentered axes - if(np.any(zerocentered==True)): - field = np.fft.fftshift(field,axes=shiftaxes(zerocentered)) - return field - - -def random_hermitian_pm1(datatype,zerocentered,shape): - - """ - Draws a set of hermitianized random, complex pm1 numbers. - - """ - - field = np.random.randint(4,high=None,size=np.prod(shape,axis=0,dtype=np.int,out=None)).reshape(shape,order='C') - dummy = np.copy(field) - ## mirror field - for ii in range(field.ndim): - dummy = np.swapaxes(dummy,0,ii) - dummy = np.flipud(dummy) - dummy = np.roll(dummy,1,axis=0) - dummy = np.swapaxes(dummy,0,ii) - field = (field+dummy+2*(field>dummy)*((field+dummy)%2))%4 ## wicked magic - x = np.array([1+0j,0+1j,-1+0j,0-1j],dtype=datatype)[field] - ## (re)shift zerocentered axes - if(np.any(zerocentered==True)): - field = np.fft.fftshift(field,axes=shiftaxes(zerocentered)) - return x - - -#----------------------------------------------------------------------------- -# Auxiliary functions -#----------------------------------------------------------------------------- - -def shiftaxes(zerocentered,st_to_zero_mode=False): - - """ - Shifts the axes in a special way needed for some functions - """ - - axes = [] - for ii in range(len(zerocentered)): - if(st_to_zero_mode==False)and(zerocentered[ii]): - axes += [ii] - if(st_to_zero_mode==True)and(not zerocentered[ii]): - axes += [ii] - return axes - - -#def nkdict(axes,dgrid,fourier=True): -# """ -# Calculates an n-dimensional array with its entries being the lengths of -# the k-vectors from the zero point of the Fourier grid. -# -# """ -# if(fourier): -# dk = dgrid -# else: -# dk = np.array([1/axes[i]/dgrid[i] for i in range(len(axes))]) -# -# kdict = np.empty(axes) -# for ii in np.ndindex(kdict.shape): -# kdict[ii] = np.sqrt(np.sum(((ii-axes//2)*dk)**2)) -# return kdict -# -# -def nkdict_fast(axes,dgrid,fourier=True): - """ - Calculates an n-dimensional array with its entries being the lengths of - the k-vectors from the zero point of the Fourier grid. - - """ - if(fourier): - dk = dgrid - else: - dk = np.array([1/dgrid[i]/axes[i] for i in range(len(axes))]) - - temp_vecs = np.array(np.where(np.ones(axes)),dtype='float').reshape(np.append(len(axes),axes)) - temp_vecs = np.rollaxis(temp_vecs,0,len(temp_vecs.shape)) - temp_vecs -= axes//2 - temp_vecs *= dk - temp_vecs *= temp_vecs - return np.sqrt(np.sum((temp_vecs),axis=-1)) - - -#def nkdict_fast2(axes,dgrid,fourier=True): -# """ -# Calculates an n-dimensional array with its entries being the lengths of -# the k-vectors from the zero point of the grid. -# -# """ -# if(fourier): -# dk = dgrid -# else: -# dk = np.array([1/dgrid[i]/axes[i] for i in range(len(axes))]) -# -# inds = [] -# for a in axes: -# inds += [slice(0,a)] -# cords = np.ogrid[inds] -# -# dists = ((cords[0]-axes[0]//2)*dk[0])**2 -# for ii in range(1,len(axes)): -# dists = dists + ((cords[ii]-axes[ii]//2)*dk[ii])**2 -# dists = np.sqrt(dists) -# -# return dists -# -# -def nklength(kdict): - return np.sort(list(set(kdict.flatten()))) - - -def drawherm(klength,kdict,ps): - - """ - Draws a hermitian random field from a Gaussian distribution. - - """ - -# vector = np.zeros(kdict.shape,dtype='complex') -# for ii in np.ndindex(vector.shape): -# if(vector[ii]==np.complex(0.,0.)): -# vector[ii] = np.sqrt(0.5*ps[np.searchsorted(klength,kdict[ii])])*np.complex(np.random.normal(0.,1.),np.random.normal(0.,1.)) -# negii = tuple(-np.array(ii)) -# vector[negii] = np.conjugate(vector[ii]) -# if(vector[negii]==vector[ii]): -# vector[ii] = np.float(np.sqrt(ps[np.searchsorted(klength,kdict[ii])]))*np.random.normal(0.,1.) -# return vector - vec = np.random.normal(loc=0,scale=1,size=kdict.size).reshape(kdict.shape) - vec = np.fft.fftn(vec)/np.sqrt(np.prod(kdict.shape)) - for ii in np.ndindex(kdict.shape): - vec[ii] *= np.sqrt(ps[np.searchsorted(klength,kdict[ii])]) - return vec - - -#def drawwild(vector,klength,kdict,ps,real_corr=1): ## vector = np.zeros(kdict.shape,dtype=np.complex) -# for ii in np.ndindex(vector.shape): -# vector[ii] = np.sqrt(real_corr*0.5*ps[klength==kdict[ii]])*np.complex(np.random.normal(0.,1.),np.random.normal(0.,1.)) -# return vector - -def drawwild(klength,kdict,ps,real_corr=1): - - """ - Draws a field of arbitrary symmetry from a Gaussian distribution. - - """ - - vec = np.empty(kdict.size,dtype=np.complex) - vec.real = np.random.normal(loc=0,scale=np.sqrt(real_corr*0.5),size=kdict.size) - vec.imag = np.random.normal(loc=0,scale=np.sqrt(real_corr*0.5),size=kdict.size) - vec = vec.reshape(kdict.shape) - for ii in np.ndindex(kdict.shape): - vec[ii] *= np.sqrt(ps[np.searchsorted(klength,kdict[ii])]) - return vec -