nifty_rg.py 95.2 KB
Newer Older
Marco Selig's avatar
Marco Selig committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2015 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/>.

"""
    ..                  __   ____   __
    ..                /__/ /   _/ /  /_
    ..      __ ___    __  /  /_  /   _/  __   __
    ..    /   _   | /  / /   _/ /  /   /  / /  /
    ..   /  / /  / /  / /  /   /  /_  /  /_/  /
    ..  /__/ /__/ /__/ /__/    \___/  \___   /  rg
    ..                               /______/

Marco Selig's avatar
Marco Selig committed
31
    NIFTY submodule for regular Cartesian grids.
Marco Selig's avatar
Marco Selig committed
32
33
34

"""
from __future__ import division
Ultimanet's avatar
Ultimanet committed
35

Marco Selig's avatar
Marco Selig committed
36
37
import os
import numpy as np
Ultimanet's avatar
Ultimanet committed
38
from scipy.special import erf 
Marco Selig's avatar
Marco Selig committed
39
40
41
import pylab as pl
from matplotlib.colors import LogNorm as ln
from matplotlib.ticker import LogFormatter as lf
Ultimanet's avatar
Ultimanet committed
42

Ultimanet's avatar
Ultimanet committed
43
44
from mpi4py import MPI

Ultimanet's avatar
Ultimanet committed
45
46
from nifty.nifty_about import about
from nifty.nifty_core import point_space,                                    \
Marco Selig's avatar
Marco Selig committed
47
                             field
Ultimanet's avatar
Ultimanet committed
48
from nifty.nifty_random import random
49
from nifty.nifty_mpi_data import distributed_data_object
Ultimanet's avatar
Ultimanet committed
50
51
52
53
from nifty.nifty_paradict import rg_space_paradict

import fft_rg

ultimanet's avatar
ultimanet committed
54
'''
Marco Selig's avatar
Marco Selig committed
55
56
57
58
59
try:
    import gfft as gf
except(ImportError):
    about.infos.cprint('INFO: "plain" gfft version 0.1.0')
    import gfft_rg as gf
ultimanet's avatar
ultimanet committed
60
'''
Ultimanet's avatar
Ultimanet committed
61
62


Marco Selig's avatar
Marco Selig committed
63
64
65
66


##-----------------------------------------------------------------------------

67
class rg_space(point_space):
Marco Selig's avatar
Marco Selig committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    """
        ..      _____   _______
        ..    /   __/ /   _   /
        ..   /  /    /  /_/  /
        ..  /__/     \____  /  space class
        ..          /______/

        NIFTY subclass for spaces of regular Cartesian grids.

        Parameters
        ----------
        num : {int, numpy.ndarray}
            Number of gridpoints or numbers of gridpoints along each axis.
        naxes : int, *optional*
            Number of axes (default: None).
        zerocenter : {bool, numpy.ndarray}, *optional*
            Whether the Fourier zero-mode is located in the center of the grid
            (or the center of each axis speparately) or not (default: True).
        hermitian : bool, *optional*
            Whether the fields living in the space follow hermitian symmetry or
            not (default: True).
        purelyreal : bool, *optional*
            Whether the field values are purely real (default: True).
        dist : {float, numpy.ndarray}, *optional*
            Distance between two grid points along each axis (default: None).
        fourier : bool, *optional*
            Whether the space represents a Fourier or a position grid
            (default: False).

        Notes
        -----
        Only even numbers of grid points per axis are supported.
        The basis transformations between position `x` and Fourier mode `k`
        rely on (inverse) fast Fourier transformations using the
        :math:`exp(2 \pi i k^\dagger x)`-formulation.

        Attributes
        ----------
        para : numpy.ndarray
            One-dimensional array containing information on the axes of the
            space in the following form: The first entries give the grid-points
            along each axis in reverse order; the next entry is 0 if the
            fields defined on the space are purely real-valued, 1 if they are
            hermitian and complex, and 2 if they are not hermitian, but
            complex-valued; the last entries hold the information on whether
            the axes are centered on zero or not, containing a one for each
            zero-centered axis and a zero for each other one, in reverse order.
        datatype : numpy.dtype
            Data type of the field values for a field defined on this space,
            either ``numpy.float64`` or ``numpy.complex128``.
        discrete : bool
            Whether or not the underlying space is discrete, always ``False``
            for regular grids.
        vol : numpy.ndarray
            One-dimensional array containing the distances between two grid
            points along each axis, in reverse order. By default, the total
            length of each axis is assumed to be one.
        fourier : bool
            Whether or not the grid represents a Fourier basis.
    """
    epsilon = 0.0001 ## relative precision for comparisons

Ultimanet's avatar
Ultimanet committed
130
    def __init__(self, num, naxes=None, zerocenter=False, hermitian=True,\
131
                purelyreal=True, dist=None, fourier=False):
Marco Selig's avatar
Marco Selig committed
132
133
134
135
136
137
138
139
140
141
142
143
        """
            Sets the attributes for an rg_space class instance.

            Parameters
            ----------
            num : {int, numpy.ndarray}
                Number of gridpoints or numbers of gridpoints along each axis.
            naxes : int, *optional*
                Number of axes (default: None).
            zerocenter : {bool, numpy.ndarray}, *optional*
                Whether the Fourier zero-mode is located in the center of the
                grid (or the center of each axis speparately) or not
Ultimanet's avatar
Ultimanet committed
144
                (default: False).
Marco Selig's avatar
Marco Selig committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
            hermitian : bool, *optional*
                Whether the fields living in the space follow hermitian
                symmetry or not (default: True).
            purelyreal : bool, *optional*
                Whether the field values are purely real (default: True).
            dist : {float, numpy.ndarray}, *optional*
                Distance between two grid points along each axis
                (default: None).
            fourier : bool, *optional*
                Whether the space represents a Fourier or a position grid
                (default: False).

            Returns
            -------
            None
        """
161

162
        complexity = 2-(bool(hermitian) or bool(purelyreal))-bool(purelyreal)
Ultimanet's avatar
Ultimanet committed
163
164
165
        if np.isscalar(num):
            num = (num,)*np.asscalar(np.array(naxes))
            
166
167
168
169
170
        self.paradict = rg_space_paradict(num=num, complexity=complexity, 
                                          zerocenter=zerocenter)        
        
        
        naxes = len(self.paradict['num'])
Marco Selig's avatar
Marco Selig committed
171
172

        ## set data type
Ultimanet's avatar
Ultimanet committed
173
        if  self.paradict['complexity'] == 0:
Marco Selig's avatar
Marco Selig committed
174
175
176
177
178
179
180
181
            self.datatype = np.float64
        else:
            self.datatype = np.complex128

        self.discrete = False

        ## set volume
        if(dist is None):
182
            dist = 1/np.array(self.paradict['num'], dtype=self.datatype)
Marco Selig's avatar
Marco Selig committed
183
        elif(np.isscalar(dist)):
184
185
            dist = self.datatype(dist)*np.ones(naxes,dtype=self.datatype,\
                                                order='C')
Marco Selig's avatar
Marco Selig committed
186
187
        else:
            dist = np.array(dist,dtype=self.datatype)
188
            if(np.size(dist) == 1):
Marco Selig's avatar
Marco Selig committed
189
190
                dist = dist*np.ones(naxes,dtype=self.datatype,order='C')
            if(np.size(dist)!=naxes):
191
192
193
                raise ValueError(about._errors.cstring(\
                    "ERROR: size mismatch ( "+str(np.size(dist))+" <> "+\
                    str(naxes)+" )."))
Marco Selig's avatar
Marco Selig committed
194
        if(np.any(dist<=0)):
195
196
197
            raise ValueError(about._errors.cstring(\
                "ERROR: nonpositive distance(s)."))
        self.vol = np.real(dist)
Marco Selig's avatar
Marco Selig committed
198
199

        self.fourier = bool(fourier)
ultimanet's avatar
ultimanet committed
200
201
        
        ## Initializes the fast-fourier-transform machine, which will be used 
ultimanet's avatar
ultimanet committed
202
        ## to transform the space
ultimanet's avatar
ultimanet committed
203
        self.fft_machine = fft_rg.fft_factory()
204
205
206
207
        
        ## Initialize the power_indices object which takes care of kindex,
        ## pindex, rho and the pundex for a given set of parameters
        if self.fourier:        
Ultimanet's avatar
Ultimanet committed
208
            self.power_indices = power_indices(shape=self.shape(),
209
210
211
                                dgrid = dist,
                                zerocentered = self.paradict['zerocenter']
                                )
Marco Selig's avatar
Marco Selig committed
212

213
214
215
216
217
218
219
220
221
222
223
224
225
    @property
    def para(self):
        temp = np.array(self.paradict['num'] + \
                         [self.paradict['complexity']] + \
                         self.paradict['zerocenter'], dtype=int)
        return temp
        
    
    @para.setter
    def para(self, x):
        self.paradict['num'] = x[:(np.size(x)-1)//2]
        self.paradict['zerocenter'] = x[(np.size(x)+1)//2:]
        self.paradict['complexity'] = x[(np.size(x)-1)//2]
Ultimanet's avatar
Ultimanet committed
226
227
228
229
230
231
232

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def apply_scalar_function(self, x, function, inplace=False):
        return x.apply_scalar_function(function, inplace=inplace)

    
233
234
235
236
237
238
239
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++      
    def unary_operation(self, x, op='None', **kwargs):
        """
        x must be a distributed_data_object which is compatible with the space!
        Valid operations are
        
        """
ultimanet's avatar
ultimanet committed
240
        
241
242
243
244
245
246
247
        translation = {"pos" : lambda y: getattr(y, '__pos__')(),
                        "neg" : lambda y: getattr(y, '__neg__')(),
                        "abs" : lambda y: getattr(y, '__abs__')(),
                        "nanmin" : lambda y: getattr(y, 'nanmin')(),
                        "min" : lambda y: getattr(y, 'amin')(),
                        "nanmax" : lambda y: getattr(y, 'nanmax')(),
                        "max" : lambda y: getattr(y, 'amax')(),
Ultimanet's avatar
Ultimanet committed
248
                        "median" : lambda y: getattr(y, 'median')(),
249
250
251
252
253
254
255
256
257
258
259
260
261
                        "mean" : lambda y: getattr(y, 'mean')(),
                        "std" : lambda y: getattr(y, 'std')(),
                        "var" : lambda y: getattr(y, 'var')(),
                        "argmin" : lambda y: getattr(y, 'argmin')(),
                        "argmin_flat" : lambda y: getattr(y, 'argmin_flat')(),
                        "argmax" : lambda y: getattr(y, 'argmax')(),
                        "argmax_flat" : lambda y: getattr(y, 'argmax_flat')(),
                        "conjugate" : lambda y: getattr(y, 'conjugate')(),
                        "None" : lambda y: y}
                        
        return translation[op](x, **kwargs)      


Marco Selig's avatar
Marco Selig committed
262
263
264
265
266
267
268
269
270
271
272
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def naxes(self):
        """
            Returns the number of axes of the grid.

            Returns
            -------
            naxes : int
                Number of axes of the regular grid.
        """
273
274
#        return (np.size(self.para)-1)//2
        return len(self.shape())
Marco Selig's avatar
Marco Selig committed
275
276
277
278
279
280
281
282
283
284

    def zerocenter(self):
        """
            Returns information on the centering of the axes.

            Returns
            -------
            zerocenter : numpy.ndarray
                Whether the grid is centered on zero for each axis or not.
        """
285
286
        #return self.para[-(np.size(self.para)-1)//2:][::-1].astype(np.bool)
        return self.paradict['zerocenter']
Marco Selig's avatar
Marco Selig committed
287
288
289
290
291
292
293
294
295
296

    def dist(self):
        """
            Returns the distances between grid points along each axis.

            Returns
            -------
            dist : np.ndarray
                Distances between two grid points on each axis.
        """
297
298
299
300
        return self.vol
 
    def shape(self):
        return np.array(self.paradict['num'])
Marco Selig's avatar
Marco Selig committed
301

Ultimanet's avatar
Ultimanet committed
302
    def dim(self, split=False):
Marco Selig's avatar
Marco Selig committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
        """
            Computes the dimension of the space, i.e.\  the number of pixels.

            Parameters
            ----------
            split : bool, *optional*
                Whether to return the dimension split up, i.e. the numbers of
                pixels along each axis, or their product (default: False).

            Returns
            -------
            dim : {int, numpy.ndarray}
                Dimension(s) of the space. If ``split==True``, a
                one-dimensional array with an entry for each axis is returned.
        """
        ## dim = product(n)
319
320
        if split == True:
            return self.shape()
Marco Selig's avatar
Marco Selig committed
321
        else:
322
            return np.prod(self.shape())
Marco Selig's avatar
Marco Selig committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def dof(self):
        """
            Computes the number of degrees of freedom of the space, i.e.\  the
            number of grid points multiplied with one or two, depending on
            complex-valuedness and hermitian symmetry of the fields.

            Returns
            -------
            dof : int
                Number of degrees of freedom of the space.
        """
        ## dof ~ dim
338
339
        if self.paradict['complexity'] < 2:
            return np.prod(self.paradict['num'])
Marco Selig's avatar
Marco Selig committed
340
        else:
341
342
343
344
345
346
            return 2*np.prod(self.paradict['num'])

#        if(self.para[(np.size(self.para)-1)//2]<2):
#            return np.prod(self.para[:(np.size(self.para)-1)//2],axis=0,dtype=None,out=None)
#        else:
#            return 2*np.prod(self.para[:(np.size(self.para)-1)//2],axis=0,dtype=None,out=None)
Marco Selig's avatar
Marco Selig committed
347

ultimanet's avatar
ultimanet committed
348

Marco Selig's avatar
Marco Selig committed
349
350
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

ultimanet's avatar
ultimanet committed
351

352
353
    def enforce_power(self, spec, size=None, kindex=None, codomain=None,
                      log=False, nbin=None, binbounds=None):
Marco Selig's avatar
Marco Selig committed
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
        """
            Provides a valid power spectrum array from a given object.

            Parameters
            ----------
            spec : {float, list, numpy.ndarray, nifty.field, function}
                Fiducial power spectrum from which a valid power spectrum is to
                be calculated. Scalars are interpreted as constant power
                spectra.

            Returns
            -------
            spec : numpy.ndarray
                Valid power spectrum.

            Other parameters
            ----------------
            size : int, *optional*
                Number of bands the power spectrum shall have (default: None).
            kindex : numpy.ndarray, *optional*
                Scale of each band.
            codomain : nifty.space, *optional*
                A compatible codomain for power indexing (default: None).
            log : bool, *optional*
378
379
380
381
                Flag specifying if the spectral binning is performed on 
                logarithmic scale or not; if set, the number of used bins is 
                set automatically (if not given otherwise); by default no 
                binning is done (default: None).
Marco Selig's avatar
Marco Selig committed
382
            nbin : integer, *optional*
383
384
385
                Number of used spectral bins; if given `log` is set to 
                ``False``; iintegers below the minimum of 3 induce an automatic
                setting; by default no binning is done (default: None).
Marco Selig's avatar
Marco Selig committed
386
387
388
            binbounds : {list, array}, *optional*
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
389
                (default: None).
Marco Selig's avatar
Marco Selig committed
390
        """
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
        
        
        
        ## Setting up the local variables: kindex 
        ## The kindex is only necessary if spec is a function or if 
        ## the size is not set explicitly 
        if kindex == None and (size == None or callable(spec) == True):
            ## Determine which space should be used to get the kindex
            if self.fourier == True:
                kindex_supply_space = self
            else:
                ## Check if the given codomain is compatible with the space  
                try:                
                    assert(self.check_codomain(codomain))
                    kindex_supply_space = codomain
                except(AssertionError):
                    about.warnings.cprint("WARNING: Supplied codomain is "+\
                    "incompatible. Generating a generic codomain. This can "+\
                    "be expensive!")
                    kindex_supply_space = self.get_codomain()
            kindex = kindex_supply_space.\
                        power_indices.get_index_dict(log=log, nbin=nbin,
                                                     binbounds=binbounds)\
                                                     ['kindex']
        
Marco Selig's avatar
Marco Selig committed
416

417
418
419
420
421
422
423
424
        
        ## Now it's about to extract a powerspectrum from spec
        ## First of all just extract a numpy array. The shape is cared about
        ## later.
                    
        ## Case 1: spec is a function
        if callable(spec) == True:
            ## Try to plug in the kindex array in the function directly            
Marco Selig's avatar
Marco Selig committed
425
            try:
426
                spec = np.array(spec(kindex), dtype=self.datatype)
Marco Selig's avatar
Marco Selig committed
427
            except:
428
429
430
431
432
433
434
435
436
437
438
439
440
441
                ## Second try: Use a vectorized version of the function.
                ## This is slower, but better than nothing
                try:
                    spec = np.vectorize(spec)(kindex)
                except:
                    raise TypeError(about._errors.cstring(
                        "ERROR: invalid power spectra function.")) 
    
        ## Case 2: spec is a field:
        elif isinstance(spec, field):
            spec = spec[:]
            spec = np.array(spec, dtype = self.datatype).flatten()
            
        ## Case 3: spec is a scalar or something else:
Marco Selig's avatar
Marco Selig committed
442
        else:
443
444
445
446
447
448
449
450
451
452
453
454
            spec = np.array(spec, dtype = self.datatype).flatten()
        
            
        ## Make some sanity checks
        ## Drop imaginary part
        temp_spec = np.real(spec)
        try:
            np.testing.assert_allclose(spec, temp_spec)
        except(AssertionError):
            about.warnings.cprint("WARNING: Dropping imaginary part.")
        spec = temp_spec
        
Marco Selig's avatar
Marco Selig committed
455
        ## check finiteness
456
        if not np.all(np.isfinite(spec)):
Marco Selig's avatar
Marco Selig committed
457
            about.warnings.cprint("WARNING: infinite value(s).")
458
        
Marco Selig's avatar
Marco Selig committed
459
        ## check positivity (excluding null)
460
461
462
463
464
465
466
        if np.any(spec<0):
            raise ValueError(about._errors.cstring(
                                "ERROR: nonpositive value(s)."))
        if np.any(spec==0):
            about.warnings.cprint("WARNING: nonpositive value(s).")            
        
        ## Set the size parameter        
Ultimanet's avatar
Ultimanet committed
467
468
        if size == None:
            size = len(kindex)
469
470
471
472
473
474
475
476
477
478
479
480
        
        ## Fix the size of the spectrum
        ## If spec is singlevalued, expand it
        if np.size(spec) == 1:
            spec = spec*np.ones(size, dtype=spec.dtype, order='C')
        ## If the size does not fit at all, throw an exception
        elif np.size(spec) < size:
            raise ValueError(about._errors.cstring("ERROR: size mismatch ( "+\
                             str(np.size(spec))+" < "+str(size)+" )."))
        elif np.size(spec) > size:
            about.warnings.cprint("WARNING: power spectrum cut to size ( == "+\
                                str(size)+" ).")
Marco Selig's avatar
Marco Selig committed
481
            spec = spec[:size]
482
        
Marco Selig's avatar
Marco Selig committed
483
484
        return spec

485

Marco Selig's avatar
Marco Selig committed
486
487
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Ultimanet's avatar
Ultimanet committed
488
    def set_power_indices(self, log=False, nbin=None, binbounds=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
        """
            Sets the (un)indexing objects for spectral indexing internally.

            Parameters
            ----------
            log : bool
                Flag specifying if the binning is performed on logarithmic
                scale or not; if set, the number of used bins is set
                automatically (if not given otherwise); by default no binning
                is done (default: None).
            nbin : integer
                Number of used bins; if given `log` is set to ``False``;
                integers below the minimum of 3 induce an automatic setting;
                by default no binning is done (default: None).
            binbounds : {list, array}
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
                (default: None).

            Returns
            -------
            None

            See Also
            --------
            get_power_indices

            Raises
            ------
            AttributeError
                If ``self.fourier == False``.
            ValueError
                If the binning leaves one or more bins empty.

        """
524
525
526
527
528

        about.warnings.cflush("WARNING: set_power_indices is a deprecated"+\
                                "function. Please use the interface of"+\
                                "self.power_indices in future!")
        self.power_indices.set_default(log=log, nbin=nbin, binbounds=binbounds)
Marco Selig's avatar
Marco Selig committed
529
530
531
        return None

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564

    def cast(self, x, verbose=False):
        """
            Computes valid field values from a given object, trying
            to translate the given data into a valid form. Thereby it is as 
            benevolent as possible. 

            Parameters
            ----------
            x : {float, numpy.ndarray, nifty.field}
                Object to be transformed into an array of valid field values.

            Returns
            -------
            x : numpy.ndarray, distributed_data_object
                Array containing the field values, which are compatible to the
                space.

            Other parameters
            ----------------
            verbose : bool, *optional*
                Whether the method should raise a warning if information is 
                lost during casting (default: False).
        """
        ## Case 1: x is a field
        if isinstance(x, field):
            if verbose:
                ## Check if the domain matches
                if(self != x.domain):
                    about.warnings.cflush(\
                    "WARNING: Getting data from foreign domain!")
            ## Extract the data, whatever it is, and cast it again
            return self.cast(x.val)
565
        
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
        ## Case 2: x is a distributed_data_object
        if isinstance(x, distributed_data_object):
            ## Check the shape
            if np.any(x.shape != self.shape()):           
                ## Check if at least the number of degrees of freedom is equal
                if x.dim() == self.dim():
                    ## If the number of dof is equal or 1, use np.reshape...
                    about.warnings.cflush(\
                    "WARNING: Trying to reshape the data. This operation is "+\
                    "expensive as it consolidates the full data!\n")
                    temp = x.get_full_data()
                    temp = np.reshape(temp, self.shape())             
                    ## ... and cast again
                    return self.cast(temp)
              
                else:
                    raise ValueError(about._errors.cstring(\
                    "ERROR: Data has incompatible shape!"))
                    
            ## Check the datatype
586
            if x.dtype < self.datatype:
587
                about.warnings.cflush(\
588
            "WARNING: Datatypes are uneqal/of conflicting precision (own: "\
Ultimanet's avatar
Ultimanet committed
589
590
591
                + str(self.datatype) + " <> foreign: " + str(x.dtype) \
                + ") and will be casted! "\
                + "Potential loss of precision!\n")
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
                temp = x.copy_empty(dtype=self.datatype)
                temp.set_local_data(x.get_local_data())
                temp.hermitian = x.hermitian
                x = temp
            
            ## Check hermitianity/reality
            if self.paradict['complexity'] == 0:
                if x.is_completely_real == False:
                    about.warnings.cflush(\
                    "WARNING: Data is not completely real. Imaginary part "+\
                    "will be discarded!\n")
                    temp = x.copy_empty()            
                    temp.set_local_data(np.real(x.get_local_data()))
                    x = temp
            
            elif self.paradict['complexity'] == 1:
                if x.hermitian == False and about.hermitianize.status:
                    about.warnings.cflush(\
                    "WARNING: Data gets hermitianized. This operation is "+\
                    "extremely expensive\n")
Ultimanet's avatar
Ultimanet committed
612
613
614
615
                    #temp = x.copy_empty()            
                    #temp.set_full_data(gp.nhermitianize_fast(x.get_full_data(), 
                    #    (False, )*len(x.shape)))
                    x = utilities.hermitianize(x)
616
617
618
619
620
621
                
            return x
                
        ## Case 3: x is something else
        ## Use general d2o casting 
        x = distributed_data_object(x, global_shape=self.shape(),\
Ultimanet's avatar
Ultimanet committed
622
            dtype=self.datatype)       
623
624
        ## Cast the d2o
        return self.cast(x)
Ultimanet's avatar
Ultimanet committed
625

Marco Selig's avatar
Marco Selig committed
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
    def enforce_values(self,x,extend=True):
        """
            Computes valid field values from a given object, taking care of
            data types, shape, and symmetry.

            Parameters
            ----------
            x : {float, numpy.ndarray, nifty.field}
                Object to be transformed into an array of valid field values.

            Returns
            -------
            x : numpy.ndarray
                Array containing the valid field values.

            Other parameters
            ----------------
            extend : bool, *optional*
                Whether a scalar is extented to a constant array or not
                (default: True).
        """
647
648
        about.warnings.cflush(\
            "WARNING: enforce_values is deprecated function. Please use self.cast")
Marco Selig's avatar
Marco Selig committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
        if(isinstance(x,field)):
            if(self==x.domain):
                if(self.datatype is not x.domain.datatype):
                    raise TypeError(about._errors.cstring("ERROR: inequal data types ( '"+str(np.result_type(self.datatype))+"' <> '"+str(np.result_type(x.domain.datatype))+"' )."))
                else:
                    x = np.copy(x.val)
            else:
                raise ValueError(about._errors.cstring("ERROR: inequal domains."))
        else:
            if(np.size(x)==1):
                if(extend):
                    x = self.datatype(x)*np.ones(self.dim(split=True),dtype=self.datatype,order='C')
                else:
                    if(np.isscalar(x)):
                        x = np.array([x],dtype=self.datatype)
                    else:
                        x = np.array(x,dtype=self.datatype)
            else:
                x = self.enforce_shape(np.array(x,dtype=self.datatype))

        ## hermitianize if ...
        if(about.hermitianize.status)and(np.size(x)!=1)and(self.para[(np.size(self.para)-1)//2]==1):
Ultimanet's avatar
Ultimanet committed
671
672
            #x = gp.nhermitianize_fast(x,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),special=False)
            x = utilities.hermitianize(x)
Marco Selig's avatar
Marco Selig committed
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
        ## check finiteness
        if(not np.all(np.isfinite(x))):
            about.warnings.cprint("WARNING: infinite value(s).")

        return x

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def get_random_values(self,**kwargs):
        """
            Generates random field values according to the specifications given
            by the parameters, taking into account possible complex-valuedness
            and hermitian symmetry.

            Returns
            -------
            x : numpy.ndarray
                Valid field values.

            Other parameters
            ----------------
            random : string, *optional*
                Specifies the probability distribution from which the random
                numbers are to be drawn.
                Supported distributions are:

                - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
                - "gau" (normal distribution with zero-mean and a given standard
                    deviation or variance)
                - "syn" (synthesizes from a given power spectrum)
                - "uni" (uniform distribution over [vmin,vmax[)

                (default: None).
            dev : float, *optional*
                Standard deviation (default: 1).
            var : float, *optional*
                Variance, overriding `dev` if both are specified
                (default: 1).
            spec : {scalar, list, numpy.ndarray, nifty.field, function}, *optional*
                Power spectrum (default: 1).
            pindex : numpy.ndarray, *optional*
                Indexing array giving the power spectrum index of each band
                (default: None).
            kindex : numpy.ndarray, *optional*
                Scale of each band (default: None).
            codomain : nifty.rg_space, *optional*
Ultimanet's avatar
Ultimanet committed
719
                A compatible codomain (default: None).
Marco Selig's avatar
Marco Selig committed
720
721
722
723
724
725
726
727
728
729
730
731
732
733
            log : bool, *optional*
                Flag specifying if the spectral binning is performed on logarithmic
                scale or not; if set, the number of used bins is set
                automatically (if not given otherwise); by default no binning
                is done (default: None).
            nbin : integer, *optional*
                Number of used spectral bins; if given `log` is set to ``False``;
                integers below the minimum of 3 induce an automatic setting;
                by default no binning is done (default: None).
            binbounds : {list, array}, *optional*
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
                (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
                Lower limit of the uniform distribution if ``random == "uni"``
Ultimanet's avatar
Ultimanet committed
734
735
                (default: 0).            
            vmin : float, *optional*
Marco Selig's avatar
Marco Selig committed
736
737
738
739
                Lower limit for a uniform distribution (default: 0).
            vmax : float, *optional*
                Upper limit for a uniform distribution (default: 1).
        """
Ultimanet's avatar
Ultimanet committed
740
741
742
743
744
745
746
747
        
        ## Parse the keyword arguments
        arg = random.parse_arguments(self,**kwargs)
        
        ## Prepare the empty distributed_data_object
        sample = distributed_data_object(global_shape=self.shape(), 
                                         dtype=self.datatype)

Ultimanet's avatar
Ultimanet committed
748
749
750
751
        ## Should the output be hermitianized?
        hermitianizeQ = about.hermitianize.status and \
                        self.paradict['complexity']          

Ultimanet's avatar
Ultimanet committed
752
        ## Case 1: uniform distribution over {-1,+1}/{1,i,-1,-i}
Ultimanet's avatar
Ultimanet committed
753
        if arg[0] == 'pm1' and hermitianizeQ == False:
Ultimanet's avatar
Ultimanet committed
754
755
756
757
            gen = lambda s: random.pm1(datatype=self.datatype,
                                       shape = s)
            sample.apply_generator(gen)
                        
Ultimanet's avatar
Ultimanet committed
758
759
760
761
762
763
764
765
766
767
768
769
770
        elif arg[0] == 'pm1' and hermitianizeQ == True:
            sample = self.get_random_values(random = 'uni', vmin=-1, vmax=1)
            local_data = sample.get_local_data()
            if issubclass(sample.dtype, np.complexfloating):
                temp_data = local_data.copy()
                local_data[temp_data.real >= 0.5] = 1
                local_data[(temp_data.real >= 0)*(temp_data.real < 0.5)] = -1
                local_data[(temp_data.real < 0)*(temp_data.imag >= 0)] = 1j
                local_data[(temp_data.real < 0)*(temp_data.imag < 0)] = -1j
            else:
                local_data[local_data >= 0] = 1
                local_data[local_data < 0] = -1
            sample.set_local_data(local_data)
Ultimanet's avatar
Ultimanet committed
771
772
773
774
775
776
            
        ## Case 2: normal distribution with zero-mean and a given standard
        ##         deviation or variance
        elif arg[0] == 'gau':
            gen = lambda s: random.gau(datatype=self.datatype,
                                       shape = s,
Ultimanet's avatar
Ultimanet committed
777
                                       mean = arg[1],
Ultimanet's avatar
Ultimanet committed
778
779
780
                                       dev = arg[2],
                                       var = arg[3])
            sample.apply_generator(gen)
Ultimanet's avatar
Ultimanet committed
781
782
783
            
            if hermitianizeQ == True:
                sample = utilities.hermitianize(sample)
Ultimanet's avatar
Ultimanet committed
784

Ultimanet's avatar
Ultimanet committed
785
786
        ## Case 3: uniform distribution
        elif arg[0] == "uni" and hermitianizeQ == False:
Ultimanet's avatar
Ultimanet committed
787
788
789
790
791
            gen = lambda s: random.uni(datatype=self.datatype,
                                       shape = s,
                                       vmin = arg[1],
                                       vmax = arg[2])
            sample.apply_generator(gen)
Ultimanet's avatar
Ultimanet committed
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
            
        elif arg[0] == "uni" and hermitianizeQ == True:
            ## For a hermitian uniform sample, generate a gaussian one
            ## and then convert it to a uniform one
            sample = self.get_random_values(random = 'gau')
            ## Use the cummulative of the gaussian, the error function in order 
            ## to transform it to a uniform distribution.
            if issubclass(sample.dtype, np.complexfloating):
                temp_func = lambda x: erf(x.real) + 1j*erf(x.imag)                  
            else:
                temp_func = lambda x: erf(x/np.sqrt(2))
            sample.apply_scalar_function(function = temp_func,
                                             inplace = True)
            
            ## Shift and stretch the uniform distribution into the given limits
            ## sample = (sample + 1)/2 * (vmax-vmin) + vmin
            vmin = arg[1]
            vmax = arg[2]            
            sample *= (vmax-vmin)/2.
            sample += 1/2.*(vmax+vmin)
            
Marco Selig's avatar
Marco Selig committed
813
814

        elif(arg[0]=="syn"):
Ultimanet's avatar
Ultimanet committed
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
            spec = arg[1]
            kpack = arg[2]
            harmonic_domain = arg[3]
            log = arg[4]
            nbin = arg[5]
            binbounds = arg[6]
            ## Check whether there is a kpack available or not.
            ## kpack is only used for computing kdict and extracting kindex
            ## If not, take kdict and kindex from the fourier_domain
            if kpack == None:
                power_indices =\
                    harmonic_domain.power_indices.get_index_dict(log = log,
                                                        nbin = nbin,
                                                        binbounds = binbounds)
                
                kindex = power_indices['kindex']
                kdict = power_indices['kdict']
                kpack = [power_indices['pindex'], power_indices['kindex']]
            else:
                kindex = kpack[1]
                kdict = harmonic_domain.power_indices.\
                    __compute_kdict_from_pindex_kindex__(kpack[0], kpack[1])           
                
Marco Selig's avatar
Marco Selig committed
838

Ultimanet's avatar
Ultimanet committed
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
            ## draw the random samples
            ## Case 1: self is a harmonic space
            if self.fourier:
                ## subcase 1: self is real
                ## -> simply generate a random field in fourier space and 
                ## weight the entries accordingly to the powerspectrum
                if self.paradict['complexity'] == 0:
                    ## set up the sample object. Overwrite the default from 
                    ## above to be sure, that the distribution strategy matches
                    ## with the one from kdict
                    sample = kdict.copy_empty(dtype = self.datatype)
                    ## set up the random number generator
                    gen = lambda s: np.random.normal(loc=0, scale=1, size=s)
                    ## apply the random number generator                    
                    sample.apply_generator(gen)
Marco Selig's avatar
Marco Selig committed
854

Ultimanet's avatar
Ultimanet committed
855
856
857
858
859
860
861
862
863
864
865
866
867
868
                
                ## subcase 2: self is hermitian but probably complex
                ## -> generate a real field (in position space) and transform
                ## it to harmonic space -> field in harmonic space is 
                ## hermitian. Now weight the modes accordingly to the 
                ## powerspectrum.
                elif self.paradict['complexity'] == 1:
                    temp_codomain = self.get_codomain()
                    ## set up the sample object. Overwrite the default from 
                    ## above to be sure, that the distribution strategy matches
                    ## with the one from kdict
                    sample = kdict.copy_empty(
                                            dtype = temp_codomain.datatype)
                    ## set up the random number generator
Ultimanet's avatar
Ultimanet committed
869

Ultimanet's avatar
Ultimanet committed
870
871
872
                    gen = lambda s: np.random.normal(loc=0, scale=1, size=s)
                    ## apply the random number generator                    
                    sample.apply_generator(gen)
Ultimanet's avatar
Ultimanet committed
873
874
875
876
877
878
879
880
881
882
                    
                    ## In order to get the normalisation right, the sqrt
                    ## of self.dim must be divided out. 
                    ## Furthermore, the normalisation in the fft routine 
                    ## must be undone
                    ## TODO: Insert explanation
                    sqrt_of_dim = np.sqrt(self.dim())
                    sample /= sqrt_of_dim
                    sample = temp_codomain.calc_weight(sample, power=-1)

Ultimanet's avatar
Ultimanet committed
883
                    ## tronsform the random field to harmonic space
Ultimanet's avatar
Ultimanet committed
884
                    sample = temp_codomain.\
Ultimanet's avatar
Ultimanet committed
885
                                        calc_transform(sample, codomain=self)
Ultimanet's avatar
Ultimanet committed
886
                    
Ultimanet's avatar
Ultimanet committed
887
888
889
890
                    ## ensure that the kdict and the harmonic_sample have the
                    ## same distribution strategy
                    assert(kdict.distribution_strategy ==\
                            sample.distribution_strategy)
Marco Selig's avatar
Marco Selig committed
891

Ultimanet's avatar
Ultimanet committed
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
                    
                ## subcase 3: self is fully complex
                ## -> generate a complex random field in harmonic space and
                ## weight the modes accordingly to the powerspectrum
                elif self.paradict['complexity'] == 2:
                    ## set up the sample object. Overwrite the default from 
                    ## above to be sure, that the distribution strategy matches
                    ## with the one from kdict
                    sample = kdict.copy_empty(dtype = self.datatype)
                    ## set up the random number generator
                    gen = lambda s: (
                        np.random.normal(loc=0, scale=1/np.sqrt(2), size=s)+
                        np.random.normal(loc=0, scale=1/np.sqrt(2), size=s)*1.j
                        )
                    ## apply the random number generator                    
                    sample.apply_generator(gen)
                
                ## apply the powerspectrum renormalization
                ## therefore extract the local data from kdict
                local_kdict = kdict.get_local_data()
                rescaler = np.sqrt(
                            spec[np.searchsorted(kindex,local_kdict)])
                sample.apply_scalar_function(lambda x: x*rescaler, 
                                             inplace=True)
            ## Case 2: self is a position space
            else:
                ## get a suitable codomain
                temp_codomain = self.get_codomain()                   

                ## subcase 1: self is a real space. 
                ## -> generate a hermitian sample with the codomain in harmonic
                ## space and make a fourier transformation.
                if self.paradict['complexity'] == 0:
                    ## check that the codomain is hermitian
                    assert(temp_codomain.paradict['complexity'] == 1)
                                                          
                ## subcase 2: self is hermitian but probably complex
                ## -> generate a real-valued random sample in fourier space
                ## and transform it to real space
                elif self.paradict['complexity'] == 1:
                    ## check that the codomain is real
                    assert(temp_codomain.paradict['complexity'] == 0)            
                
                ## subcase 3: self is fully complex
                ## -> generate a complex-valued random sample in fourier space
                ## and transform it to real space
                elif self.paradict['complexity'] == 2:
                    ## check that the codomain is real
                    assert(temp_codomain.paradict['complexity'] == 2)            


                ## Get a hermitian/real/complex sample in harmonic space from 
                ## the codomain
                sample = temp_codomain.get_random_values(
                                                    random='syn',
                                                    pindex = kpack[0],
                                                    kindex = kpack[1],
                                                    spec = spec,
                                                    codomain = self,
                                                    log = log,
                                                    nbin = nbin,
                                                    binbounds = binbounds
                                                    )
Ultimanet's avatar
Ultimanet committed
955
956
957
                ## Correct the weighting
                #sample = self.calc_weight(sample, power=-1)
                
Ultimanet's avatar
Ultimanet committed
958
959
                ## Take the fourier transform
                sample = temp_codomain.calc_transform(sample, 
Ultimanet's avatar
Ultimanet committed
960
                                                      codomain = self)
Ultimanet's avatar
Ultimanet committed
961
962
963
964
965
966
967
968

            if self.paradict['complexity'] == 1:
               sample.hermitian = True
           
        else:
            raise KeyError(about._errors.cstring(
                        "ERROR: unsupported random key '"+str(arg[0])+"'."))
     
Ultimanet's avatar
Ultimanet committed
969
             
Ultimanet's avatar
Ultimanet committed
970
971
972
973
        return sample



Marco Selig's avatar
Marco Selig committed
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991


    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def check_codomain(self,codomain):
        """
            Checks whether a given codomain is compatible to the space or not.

            Parameters
            ----------
            codomain : nifty.space
                Space to be checked for compatibility.

            Returns
            -------
            check : bool
                Whether or not the given codomain is compatible to the space.
        """
992
993
994
        if codomain == None:
            return False
            
Ultimanet's avatar
Ultimanet committed
995
        if(not isinstance(codomain,rg_space)):
Marco Selig's avatar
Marco Selig committed
996
997
            raise TypeError(about._errors.cstring("ERROR: invalid input."))

Ultimanet's avatar
Ultimanet committed
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
        ## check number of number and size of axes 
        if not np.all(self.paradict['num'] == codomain.paradict['num']):
            return False
            
        ## check fourier flag
        if self.fourier == codomain.fourier:
            return False
            
        ## check complexity-type
        ## prepare the shorthands
        dcomp = self.paradict['complexity']
        cocomp = codomain.paradict['complexity']
        
        ## Case 1: if the domain is copmleteley complex 
        ## -> the codomain must be complex, too
        if dcomp == 2:
            if cocomp != 2:
                return False
        ## Case 2: domain is hermitian
        ## -> codmomain can be real. If it is marked as hermitian or even 
        ## fully complex, a warning is raised
        elif dcomp == 1:
            if cocomp > 0:
                about.warnings.cprint("WARNING: Unrecommended codomain! "+\
                    "The domain is hermitian, hence the codomain should "+\
                    "be restricted to real values!")
        
        ## Case 3: domain is real
        ## -> codmain should be hermitian
        elif dcomp == 0:
            if cocomp == 2:
                about.warnings.cprint("WARNING: Unrecommended codomain! "+\
                    "The domain is real, hence the codomain should "+\
                    "be restricted to hermitian configurations!")
            elif cocomp == 0:
                return False

        ## Check if the distances match, i.e. dist'=1/(num*dist)
        if not np.all(
                np.absolute(self.paradict['num']*
                            self.vol*
                            codomain.vol-1) < self.epsilon):
            return False
            
        return True
Marco Selig's avatar
Marco Selig committed
1043

Ultimanet's avatar
Ultimanet committed
1044
        
Marco Selig's avatar
Marco Selig committed
1045
1046
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Ultimanet's avatar
Ultimanet committed
1047
    def get_codomain(self, coname=None, cozerocenter=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
        """
            Generates a compatible codomain to which transformations are
            reasonable, i.e.\  either a shifted grid or a Fourier conjugate
            grid.

            Parameters
            ----------
            coname : string, *optional*
                String specifying a desired codomain (default: None).
            cozerocenter : {bool, numpy.ndarray}, *optional*
                Whether or not the grid is zerocentered for each axis or not
                (default: None).

            Returns
            -------
            codomain : nifty.rg_space
                A compatible codomain.

            Notes
            -----
            Possible arguments for `coname` are ``'f'`` in which case the
Ultimanet's avatar
Ultimanet committed
1069
1070
1071
            codomain arises from a Fourier transformation, ``'i'`` in which 
            case it arises from an inverse Fourier transformation.If no 
            `coname` is given, the Fourier conjugate grid is produced.
Marco Selig's avatar
Marco Selig committed
1072
        """
Ultimanet's avatar
Ultimanet committed
1073
1074
        naxes = self.naxes() 
        ## Parse the cozerocenter input
Marco Selig's avatar
Marco Selig committed
1075
        if(cozerocenter is None):
Ultimanet's avatar
Ultimanet committed
1076
1077
            cozerocenter = self.paradict['zerocenter']
        ## if the input is something scalar, cast it to a boolean
Marco Selig's avatar
Marco Selig committed
1078
1079
        elif(np.isscalar(cozerocenter)):
            cozerocenter = bool(cozerocenter)
Ultimanet's avatar
Ultimanet committed
1080
        ## if it is not a scalar...
Marco Selig's avatar
Marco Selig committed
1081
        else:
Ultimanet's avatar
Ultimanet committed
1082
            ## ...cast it to a numpy array of booleans
Marco Selig's avatar
Marco Selig committed
1083
            cozerocenter = np.array(cozerocenter,dtype=np.bool)
Ultimanet's avatar
Ultimanet committed
1084
            ## if it was a list of length 1, extract the boolean
Marco Selig's avatar
Marco Selig committed
1085
1086
            if(np.size(cozerocenter)==1):
                cozerocenter = np.asscalar(cozerocenter)
Ultimanet's avatar
Ultimanet committed
1087
1088
            ## if the length of the input does not match the number of 
            ## dimensions, raise an exception
Marco Selig's avatar
Marco Selig committed
1089
            elif(np.size(cozerocenter)!=naxes):
Ultimanet's avatar
Ultimanet committed
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
                raise ValueError(about._errors.cstring(
                    "ERROR: size mismatch ( "+\
                    str(np.size(cozerocenter))+" <> "+str(naxes)+" )."))
        
        ## Set up the initialization variables
        num = self.paradict['num']
        purelyreal = (self.paradict['complexity'] == 1)        
        hermitian = (self.paradict['complexity'] < 2)
        dist = 1/(self.paradict['num']*self.vol)        
        
        if coname == None:
            fourier = bool(not self.fourier)            
        elif coname[0] == 'f':
            fourier = True
        elif coname[0] == 'i':
            fourier = False
        else:
            raise ValueError(about._errors.cstring(
                                            "ERROR: Unknown coname keyword"))
        new_space = rg_space(num,
                             zerocenter = cozerocenter,
                             hermitian = hermitian,
                             purelyreal = purelyreal,
                             dist = dist,
                             fourier = fourier)                                            
        return new_space
Marco Selig's avatar
Marco Selig committed
1116
1117
1118

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Ultimanet's avatar
Ultimanet committed
1119
    def get_meta_volume(self, total=False):
Marco Selig's avatar
Marco Selig committed
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
        """
            Calculates the meta volumes.

            The meta volumes are the volumes associated with each component of
            a field, taking into account field components that are not
            explicitly included in the array of field values but are determined
            by symmetry conditions. In the case of an :py:class:`rg_space`, the
            meta volumes are simply the pixel volumes.

            Parameters
            ----------
            total : bool, *optional*
                Whether to return the total meta volume of the space or the
                individual ones of each pixel (default: False).

            Returns
            -------
            mol : {numpy.ndarray, float}
                Meta volume of the pixels or the complete space.
        """
        if(total):
Ultimanet's avatar
Ultimanet committed
1141
            return self.dim(split=False)*np.prod(self.vol)
Marco Selig's avatar
Marco Selig committed
1142
        else:
Ultimanet's avatar
Ultimanet committed
1143
            mol = np.ones(self.dim(split=True),dtype=self.vol.dtype)
Marco Selig's avatar
Marco Selig committed
1144
1145
1146
1147
            return self.calc_weight(mol,power=1)

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Ultimanet's avatar
Ultimanet committed
1148
    def calc_weight(self, x, power=1):
Marco Selig's avatar
Marco Selig committed
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
        """
            Weights a given array with the pixel volumes to a given power.

            Parameters
            ----------
            x : numpy.ndarray
                Array to be weighted.
            power : float, *optional*
                Power of the pixel volumes to be used (default: 1).

            Returns
            -------
            y : numpy.ndarray
                Weighted array.
        """
1164
1165
1166
        #x = self.cast(x)
        if isinstance(x, distributed_data_object):
            is_hermitianQ = x.hermitian
Marco Selig's avatar
Marco Selig committed
1167
        ## weight
1168
1169
1170
        x =  x * self.get_weight(power = power)
        if isinstance(x, distributed_data_object):        
            x.hermitian = is_hermitianQ
Ultimanet's avatar
Ultimanet committed
1171
        return x
Marco Selig's avatar
Marco Selig committed
1172

1173
1174
1175
    def get_weight(self, power = 1):
        return np.prod(self.vol)**power

Marco Selig's avatar
Marco Selig committed
1176
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1177
    def calc_dot(self, x, y):
Marco Selig's avatar
Marco Selig committed
1178
        """
1179
1180
            Computes the discrete inner product of two given arrays of field
            values.
Marco Selig's avatar
Marco Selig committed
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193

            Parameters
            ----------
            x : numpy.ndarray
                First array
            y : numpy.ndarray
                Second array

            Returns
            -------
            dot : scalar
                Inner product of the two arrays.
        """
1194
1195
1196
1197
1198
        x = self.cast(x)
        y = self.cast(y)
        result = x.vdot(y)
        if np.isreal(result):
            result = np.asscalar(np.real(result))      
Ultimanet's avatar
Ultimanet committed
1199
        if self.paradict['complexity'] != 2:
1200
1201
            if(np.absolute(result.imag) > self.epsilon**2\
                                          *np.absolute(result.real)):
Ultimanet's avatar
Ultimanet committed
1202
1203
                about.warnings.cprint(
                    "WARNING: Discarding considerable imaginary part.")
1204
1205
1206
            result = np.asscalar(np.real(result))      
        return result
        
Marco Selig's avatar
Marco Selig committed
1207
1208
1209

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

1210
    def calc_transform(self, x, codomain=None, **kwargs):
Marco Selig's avatar
Marco Selig committed
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
        """
            Computes the transform of a given array of field values.

            Parameters
            ----------
            x : numpy.ndarray
                Array to be transformed.
            codomain : nifty.rg_space, *optional*
                Target space to which the transformation shall map
                (default: None).

            Returns
            -------
            Tx : numpy.ndarray
                Transformed array
        """
1227
1228
1229
        x = self.cast(x)
        
        if(codomain is None):
Ultimanet's avatar
Ultimanet committed
1230
            codomain = self.get_codomain()
1231
1232
1233
1234
1235
1236
1237
1238
1239
        
        ## Check if the given codomain is suitable for the transformation
        if (not isinstance(codomain, rg_space)) or \
                (not self.check_codomain(codomain)):
            raise ValueError(about._errors.cstring(
                                "ERROR: unsupported codomain."))
        if codomain.fourier == True:
            ## correct for forward fft
            x = self.calc_weight(x, power=1)
Ultimanet's avatar
Ultimanet committed
1240
1241
1242
1243
        
#            ## correct for inverse fft
#            x = self.calc_weight(x, power=1)
#            x *= self.dim(split=False)
1244
1245
1246
1247
1248
        
        ## Perform the transformation
        Tx = self.fft_machine.transform(val=x, domain=self, codomain=codomain, 
                                        **kwargs)

Ultimanet's avatar
Ultimanet committed
1249
1250
1251
1252
1253
        if codomain.fourier == False:
            ## correct for inverse fft
            Tx = codomain.calc_weight(Tx, power=-1)


1254
1255
1256
        ## when the target space is purely real, the result of the 
        ## transformation must be corrected accordingly. Using the casting 
        ## method of codomain is sufficient
Ultimanet's avatar
Ultimanet committed
1257
        ## TODO: Let .transform  yield the correct datatype
1258
1259
1260
1261
        Tx = codomain.cast(Tx)
        
        return Tx

Marco Selig's avatar
Marco Selig committed
1262
1263
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Ultimanet's avatar
Ultimanet committed
1264
    def calc_smooth(self, x, sigma=0, codomain=None):
Marco Selig's avatar
Marco Selig committed
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
        """
            Smoothes an array of field values by convolution with a Gaussian
            kernel.

            Parameters
            ----------
            x : numpy.ndarray
                Array of field values to be smoothed.
            sigma : float, *optional*
                Standard deviation of the Gaussian kernel, specified in units
                of length in position space; for testing: a sigma of -1 will be
                reset to a reasonable value (default: 0).

            Returns
            -------
            Gx : numpy.ndarray
                Smoothed array.
        """

Ultimanet's avatar
Ultimanet committed
1284
1285
1286
        
        ## Check sigma
        if sigma == 0:
Marco Selig's avatar
Marco Selig committed
1287
            return x
Ultimanet's avatar
Ultimanet committed
1288
1289
1290
1291
        elif sigma == -1:
            about.infos.cprint(
                "INFO: Resetting sigma to sqrt(2)*max(dist).")
            sigma = np.sqrt(2)*np.max(self.dist())
Marco Selig's avatar
Marco Selig committed
1292
1293
        elif(sigma<0):
            raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
Ultimanet's avatar
Ultimanet committed
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
        

        ## if a codomain was given...
        if codomain != None:
            ## ...check if it was suitable
            if not isinstance(codomain, rg_space):
                raise ValueError(about._errors.cstring(
                    "ERROR: codomain is not a rg_space instance!"))
            if self.fourier == False and codomain.fourier == False:
                raise ValueError(about._errors.cstring(
                    "ERROR: fourier_domain is not a fourier space!"))
            if not self.check_codomain(codomain):
                raise ValueError(about._errors.cstring(
                    "ERROR: fourier_codomain is not a valid codomain!"))
        elif self.fourier == False:
            codomain = self.get_codomain()

        ## Case1: 
        ## If self is a position-space, fourier transform the input and
        ## call calc_smooth of the fourier codomain
        if self.fourier == False:
            x = self.calc_transform(x, codomain = codomain)
            x = codomain.calc_smooth(x, sigma)
            x = codomain.calc_transform(x, codomain = self)
            return x
        
        ## Case 2: 
        ## if self is fourier multiply the gaussian kernel, etc...
        
        ## Cast the input
        x = self.cast(x)       
         
        ## if x is hermitian it remains hermitian during smoothing
        remeber_hermitianQ = x.hermitian

        ## Define the Gaussian kernel function 
        gaussian = lambda x: np.exp(-2.*np.pi**2*x**2*sigma**2)
        
        ## Define the variables in the dialect of the legacy smoothing.py 
        nx = self.shape()
        dx = 1/nx/self.vol
        ## Multiply the data along each axis with suitable the gaussian kernel
        for i in range(len(nx)):
            ## Prepare the exponent
            dk = 1./nx[i]/dx[i]
            nk = nx[i]
            k = -0.5*nk*dk + np.arange(nk)*dk
            if self.paradict['zerocenter'][i] == False:
                k = np.fft.fftshift(k)
            ## compute the actual kernel vector
            gaussian_kernel_vector = gaussian(k)
            ## blow up the vector to an array of shape (1,.,1,len(nk),1,.,1)
            blown_up_shape = [1,]*len(nx)
            blown_up_shape[i] = len(gaussian_kernel_vector)
            gaussian_kernel_vector =\
                gaussian_kernel_vector.reshape(blown_up_shape)
            ## apply the blown-up gaussian_kernel_vector
            x *= gaussian_kernel_vector
        x.hermitian = remeber_hermitianQ
        
        return x
        
Marco Selig's avatar
Marco Selig committed
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399

    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    def calc_power(self,x,**kwargs):
        """
            Computes the power of an array of field values.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values of which the power is to be
                calculated.

            Returns
            -------
            spec : numpy.ndarray
                Power contained in the input array.

            Other parameters
            ----------------
            pindex : numpy.ndarray, *optional*
                Indexing array assigning the input array components to
                components of the power spectrum (default: None).
            rho : numpy.ndarray, *optional*
                Number of degrees of freedom per band (default: None).
            codomain : nifty.space, *optional*
                A compatible codomain for power indexing (default: None).
            log : bool, *optional*
                Flag specifying if the spectral binning is performed on logarithmic
                scale or not; if set, the number of used bins is set
                automatically (if not given otherwise); by default no binning
                is done (default: None).
            nbin : integer, *optional*
                Number of used spectral bins; if given `log` is set to ``False``;
                integers below the minimum of 3 induce an automatic setting;
                by default no binning is done (default: None).
            binbounds : {list, array}, *optional*
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
                (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
                Lower limit of the uniform distribution if ``random == "uni"``
                (default: 0).

        """
Ultimanet's avatar
Ultimanet committed
1400
1401
1402
1403
        x = self.cast(x)

        ## If self is a position space, delegate calc_power to its codomain.
        if self.fourier == False:
Marco Selig's avatar
Marco Selig committed
1404
            try:
Ultimanet's avatar
Ultimanet committed
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
                codomain = kwargs.get('codomain')
            except(KeyError):
                codomain = self.get_codomain()
                
            y = self.calc_transform(x, codomain)
            kwargs.update({'codomain': self})
            return codomain.calc_power(y, **kwargs)
            
        ## If some of the pindex, kindex or rho arrays are given explicitly,
        ## favor them over those from the self.power_indices dictionary.
        ## As the default value in kwargs.get(key, default) does NOT evaluate
        ## lazy, a distinction of cases is necessary. Otherwise the 
        ## powerindices might be computed, although not necessary
        if kwargs.has_key('pindex') and kwargs.has_key('kindex') and\
                kwargs.has_key('rho'):
            pindex = kwargs.get('pindex')
            rho = kwargs.get('rho')
        else:
            log = kwargs.get('log', None)
            nbin = kwargs.get('nbin', None)
            binbounds = kwargs.get('binbounds', None)            
            power_indices = self.power_indices.get_index_dict(log = log,
                                                              nbin = nbin,
                                                        binbounds = binbounds)            
            pindex = kwargs.get('pindex', power_indices['pindex'])
            rho = kwargs.get('rho', power_indices['rho'])
        
        fieldabs = abs(x)**2
        power_spectrum = np.zeros(rho.shape) 
1434

Ultimanet's avatar
Ultimanet committed
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
                
        ## In order to make the summation over identical pindices fast, 
        ## the pindex and the kindex must have the same distribution strategy
        if pindex.distribution_strategy == fieldabs.distribution_strategy and\
            pindex.distributor.comm == fieldabs.distributor.comm:
            working_field = fieldabs
        else:
            working_field = pindex.copy_empty(dtype = fieldabs.dtype)
            working_field.inject((slice(None),), fieldabs, (slice(None,)))
        
        local_power_spectrum = np.bincount(pindex.get_local_data().flatten(),
                        weights = working_field.get_local_data().flatten())        
        power_spectrum =\
            pindex.distributor._allgather(local_power_spectrum)
        power_spectrum = np.sum(power_spectrum, axis = 0)
1450

Ultimanet's avatar
Ultimanet committed
1451
1452
1453
1454
        ## Divide out the degeneracy factor        
        power_spectrum /= rho
        return power_spectrum
        
Marco Selig's avatar
Marco Selig committed
1455

Ultimanet's avatar
Ultimanet committed
1456
        
Marco Selig's avatar
Marco Selig committed
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
    ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    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.

            Parameters
            ----------
            x : numpy.ndarray
                Array containing the field values.

            Returns
            -------
            None

            Other parameters
            ----------------
            title : string, *optional*
                Title of the plot (default: "").
            vmin : float, *optional*
                Minimum value to be displayed (default: ``min(x)``).
            vmax : float, *optional*
                Maximum value to be displayed (default: ``max(x)``).
            power : bool, *optional*
                Whether to plot the power contained in the field or the field
                values themselves (default: False).
            unit : string, *optional*
                Unit of the field values (default: "").
            norm : string, *optional*
                Scaling of the field values before plotting (default: None).
            cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
                Color map to be used for two-dimensional plots (default: None).
            cbar : bool, *optional*
                Whether to show the color bar or not (default: True).
            other : {single object, tuple of objects}, *optional*
                Object or tuple of objects to be added, where objects can be
                scalars, arrays, or fields (default: None).
            legend : bool, *optional*
                Whether to show the legend or not (default: False).
            mono : bool, *optional*
                Whether to plot the monopole or not (default: True).
            save : string, *optional*
                Valid file name where the figure is to be stored, by default
                the figure is not saved (default: False).
            error : {float, numpy.ndarray, nifty.field}, *optional*
                Object indicating some confidence interval to be plotted
                (default: None).
            kindex : numpy.ndarray, *optional*
                Scale corresponding to each band in the power spectrum
                (default: None).
            codomain : nifty.space, *optional*
                A compatible codomain for power indexing (default: None).
            log : bool, *optional*
                Flag specifying if the spectral binning is performed on logarithmic
                scale or not; if set, the number of used bins is set
                automatically (if not given otherwise); by default no binning
                is done (default: None).
            nbin : integer, *optional*
                Number of used spectral bins; if given `log` is set to ``False``;
                integers below the minimum of 3 induce an automatic setting;
                by default no binning is done (default: None).
            binbounds : {list, array}, *optional*
                User specific inner boundaries of the bins, which are preferred
                over the above parameters; by default no binning is done
                (default: None).            vmin : {scalar, list, ndarray, field}, *optional*
                Lower limit of the uniform distribution if ``random == "uni"``
                (default: 0).

        """
        if(not pl.isinteractive())and(not bool(kwargs.get("save",False))):
            about.warnings.cprint("WARNING: interactive mode off.")

        naxes = (np.size(self.para)-1)//2
        if(power is None):
            power = bool(self.para[naxes])

        if(power):
            x = self.calc_power(x,**kwargs)

            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
            xaxes = kwargs.get("kindex",None)
            ## implicit kindex
            if(xaxes is None):
                try:
                    self.set_power_indices(**kwargs)
                except:
                    codomain = kwargs.get("codomain",self.get_codomain())
                    codomain.set_power_indices(**kwargs)
                    xaxes = codomain.power_indices.get("kindex")
                else:
                    xaxes = self.power_indices.get("kindex")

            if(norm is None)or(not isinstance(norm,int)):
                norm = naxes
            if(vmin is None):
                vmin = np.min(x[:mono].tolist()+(xaxes**norm*x)[1:].tolist(),axis=None,out=None)
            if(vmax is None):
                vmax = np.max(x[:mono].tolist()+(xaxes**norm*x)[1:].tolist(),axis=None,out=None)
            ax0.loglog(xaxes[1:],(xaxes**norm*x)[1:],color=[0.0,0.5,0.0],label="graph 0",linestyle='-',linewidth=2.0,zorder=1)
            if(mono):
                ax0.scatter(0.5*(xaxes[1]+xaxes[2]),x[0],s=20,color=[0.0,0.5,0.0],marker='o',cmap=None,norm=None,vmin=None,vmax=None,alpha=None,linewidths=None,verts=None,zorder=1)

            if(other is not None):
                if(isinstance(other,tuple)):
                    other = list(other)
                    for ii in xrange(len(other)):
                        if(isinstance(other[ii],field)):
                            other[ii] = other[ii].power(**kwargs)
                        else:
                            other[ii] = self.enforce_power(other[ii],size=np.size(xaxes),kindex=xaxes)
                elif(isinstance(other,field)):
                    other = [other.power(**kwargs)]
                else:
                    other = [self.enforce_power(other,size=np.size(xaxes),kindex=xaxes)]
                imax = max(1,len(other)-1)
                for ii in xrange(len(other)):
                    ax0.loglog(xaxes[1:],(xaxes**norm*other[ii])[1:],color=[max(0.0,1.0-(2*ii/imax)**2),0.5*((2*ii-imax)/imax)**2,max(0.0,1.0-(2*(ii-imax)/imax)**2)],label="graph "+str(ii+1),linestyle='-',linewidth=1.0,zorder=-ii)
                    if(mono):
                        ax0.scatter(0.5*(xaxes[1]+xaxes[2]),other[ii][0],s=20,color=[max(0.0,1.0-(2*ii/imax)**2),0.5*((2*ii-imax)/imax)**2,max(0.0,1.0-(2*(ii-imax)/imax)**2)],marker='o',cmap=None,norm=None,vmin=None,vmax=None,alpha=None,linewidths=None,verts=None,zorder=-ii)
                if(legend):
                    ax0.legend()

            ax0.set_xlim(xaxes[1],xaxes[-1])
            ax0.set_xlabel(r"$|k|$")
            ax0.set_ylim(vmin,vmax)
            ax0.set_ylabel(r"$|k|^{%i} P_k$"%norm)
            ax0.set_title(title)

        else:
            x = self.enforce_shape(np.array(x))

            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])

                xaxes = (np.arange(self.para[0],dtype=np.int)+self.para[2]*(self.para[0]//2))*self.vol
                if(vmin is None):
                    if(np.iscomplexobj(x)):
                        vmin = min(np.min(np.absolute(x),axis=None,out=None),np.min(np.real(x),axis=None,out=None),np.min(np.imag(x),axis=None,out=None))
                    else:
                        vmin = np.min(x,axis=None,out=None)
                if(vmax is None):
                    if(np.iscomplexobj(x)):
                        vmax = max(np.max(np.absolute(x),axis=None,out=None),np.max(np.real(x),axis=None,out=None),np.max(np.imag(x),axis=None,out=None))
                    else:
                        vmax = np.max(x,axis=None,out=None)
                if(norm=="log"):
                    ax0graph = ax0.semilogy
                    if(vmin<=0):
                        raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
                else:
                    ax0graph = ax0.plot

                if(np.iscomplexobj(x)):
                    ax0graph(xaxes,np.absolute(x),color=[0.0,0.5,0.0],label="graph (absolute)",linestyle='-',linewidth=2.0,zorder=1)
                    ax0graph(xaxes,np.real(x),color=[0.0,0.5,0.0],label="graph (real part)",linestyle="--",linewidth=1.0,zorder=0)
                    ax0graph(xaxes,np.imag(x),color=[0.0,0.5,0.0],label="graph (imaginary part)",linestyle=':',linewidth=1.0,zorder=0)
                    if(legend):
                        ax0.legend()
                elif(other is not None):
                    ax0graph(xaxes,x,color=[0.0,0.5,0.0],label="graph 0",linestyle='-',linewidth=2.0,zorder=1)
                    if(isinstance(other,tuple)):
                        other = [self.enforce_values(xx,extend=True) for xx in other]
                    else:
                        other = [self.enforce_values(other,extend=True)]
                    imax = max(1,len(other)-1)
                    for ii in xrange(len(other)):
                        ax0graph(xaxes,other[ii],color=[max(0.0,1.0-(2*ii/imax)**2),0.5*((2*ii-imax)/imax)**2,max(0.0,1.0-(2*(ii-imax)/imax)**2)],label="graph "+str(ii+1),linestyle='-',linewidth=1.0,zorder=-ii)
                    if("error" in kwargs):
                        error = self.enforce_values(np.absolute(kwargs.get("error")),extend=True)
                        ax0.fill_between(xaxes,x-error,x+error,color=[0.8,0.8,0.8],label="error 0",zorder=-len(other))
                    if(legend):
                        ax0.legend()
                else:
                    ax0graph(xaxes,x,color=[0.0,0.5,0.0],label="graph 0",linestyle='-',linewidth=2.0,zorder=1)
                    if("error" in kwargs):
                        error = self.enforce_values(np.absolute(kwargs.get("error")),extend=True)
                        ax0.fill_between(xaxes,x-error,x+error,color=[0.8,0.8,0.8],label="error 0",zorder=0)

                ax0.set_xlim(xaxes[0],xaxes[-1])
                ax0.set_xlabel("coordinate")
                ax0.set_ylim(vmin,vmax)
                if(unit):
                    unit = " ["+unit+"]"
                ax0.set_ylabel("values"+unit)
                ax0.set_title(title)

            elif(naxes==2):
                if(np.