nifty_probing.py 19.6 KB
Newer Older
Ultima's avatar
Ultima committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# NIFTY (Numerical Information Field Theory) has been developed at the
# Max-Planck-Institute for Astrophysics.
#
# Copyright (C) 2015 Max-Planck-Society
#
# Author: Theo Steininger
# Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
Ultimanet's avatar
Ultimanet committed
21
22
23

from __future__ import division

24
from nifty.config import about
csongor's avatar
csongor committed
25
26
from nifty.nifty_core import space
from nifty.nifty_field import field
27
from nifty.nifty_utilities import direct_vdot
Ultimanet's avatar
Ultimanet committed
28
29


Ultima's avatar
Ultima committed
30

Ultimanet's avatar
Ultimanet committed
31

Ultima's avatar
Ultima committed
32
class prober(object):
Ultimanet's avatar
Ultimanet committed
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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
    """
        ..                                    __        __
        ..                                  /  /      /__/
        ..      ______    _____   ______   /  /___    __   __ ___    ____ __
        ..    /   _   | /   __/ /   _   | /   _   | /  / /   _   | /   _   /
        ..   /  /_/  / /  /    /  /_/  / /  /_/  / /  / /  / /  / /  /_/  /
        ..  /   ____/ /__/     \______/  \______/ /__/ /__/ /__/  \____  /  class
        .. /__/                                                  /______/

        NIFTY class for probing (using multiprocessing)

        This is the base NIFTY probing class from which other probing classes
        (e.g. diagonal probing) are derived.

        When called, a probing class instance evaluates an operator or a
        function using random fields, whose components are random variables
        with mean 0 and variance 1. When an instance is called it returns the
        mean value of f(probe), where probe is a random field with mean 0 and
        variance 1. The mean is calculated as 1/N Sum[ f(probe_i) ].

        Parameters
        ----------
        op : operator
            The operator specified by `op` is the operator to be probed.
            If no operator is given, then probing will be done by applying
            `function` to the probes. (default: None)
        function : function, *optional*
            If no operator has been specified as `op`, then specification of
            `function` is non optional. This is the function, that is applied
            to the probes. (default: `op.times`)
        domain : space, *optional*
            If no operator has been specified as `op`, then specification of
            `domain` is non optional. This is the space that the probes live
            in. (default: `op.domain`)
        target : domain, *optional*
            `target` is the codomain of `domain`
            (default: `op.domain.get_codomain()`)
        random : string, *optional*
            the distribution from which the probes are drawn. `random` can be
            either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
            or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
            zero-mean and unit-variance (default: "pm1")
        ncpu : int, *optional*
            the number of cpus to be used from parallel probing. (default: 2)
        nrun : int, *optional*
            the number of probes to be evaluated. If `nrun<ncpu`, it will be
            set to `ncpu`. (default: 8)
        nper : int, *optional*
            this number specifies how many probes will be evaluated by one
            worker. Afterwards a new worker will be created to evaluate a chunk
            of `nper` probes.
            If for example `nper=nrun/ncpu`, then every worker will be created
            for every cpu. This can lead to the case, that all workers but one
            are already finished, but have to wait for the last worker that
            might still have a considerable amount of evaluations left. This is
            obviously not very effective.
            If on the other hand `nper=1`, then for each evaluation a worker will
            be created. In this case all cpus will work until nrun probes have
            been evaluated.
            It is recommended to leave `nper` as the default value. (default: 1)
        var : bool, *optional*
            If `var` is True, then the variance of the sampled function will
            also be returned. The result is then a tuple with the mean in the
            zeroth entry and the variance in the first entry. (default: False)


        See Also
        --------
        diagonal_probing : A probing class to get the diagonal of an operator
        trace_probing : A probing class to get the trace of an operator


        Attributes
        ----------
        function : function
            the function, that is applied to the probes
        domain : space
            the space, where the probes live in
        target : space
            the codomain of `domain`
        random : string
            the random number generator used to create the probes
            (either "pm1" or "gau")
        ncpu : int
            the number of cpus used for probing
        nrun : int
            the number of probes to be evaluated, when the instance is called
        nper : int
            number of probes, that will be evaluated by one worker
        var : bool
            whether the variance will be additionally returned, when the
            instance is called
        quargs : dict
            Keyword arguments passed to `function` in each call.

    """
129
    def __init__(self, operator=None, function=None, domain=None,
Ultima's avatar
Ultima committed
130
                 codomain=None, random="pm1", nrun=8, varQ=False, **kwargs):
Ultimanet's avatar
Ultimanet committed
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
        """
        initializes a probing instance

        Parameters
        ----------
        op : operator
            The operator specified by `op` is the operator to be probed.
            If no operator is given, then probing will be done by applying
            `function` to the probes. (default: None)
        function : function, *optional*
            If no operator has been specified as `op`, then specification of
            `function` is non optional. This is the function, that is applied
            to the probes. (default: `op.times`)
        domain : space, *optional*
            If no operator has been specified as `op`, then specification of
            `domain` is non optional. This is the space that the probes live
            in. (default: `op.domain`)
        target : domain, *optional*
            `target` is the codomain of `domain`
            (default: `op.domain.get_codomain()`)
        random : string, *optional*
            the distribution from which the probes are drawn. `random` can be
            either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
            or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
            zero-mean and unit-variance (default: "pm1")
        ncpu : int, *optional*
            the number of cpus to be used from parallel probing. (default: 2)
        nrun : int, *optional*
            the number of probes to be evaluated. If `nrun<ncpu`, it will be
            set to `ncpu`. (default: 8)
        nper : int, *optional*
            this number specifies how many probes will be evaluated by one
            worker. Afterwards a new worker will be created to evaluate a chunk
            of `nper` probes.
            If for example `nper=nrun/ncpu`, then every worker will be created
            for every cpu. This can lead to the case, that all workers but one
            are already finished, but have to wait for the last worker that
            might still have a considerable amount of evaluations left. This is
            obviously not very effective.
            If on the other hand `nper=1`, then for each evaluation a worker will
            be created. In this case all cpus will work until nrun probes have
            been evaluated.
            It is recommended to leave `nper` as the default value. (default: 1)
        var : bool, *optional*
            If `var` is True, then the variance of the sampled function will
            also be returned. The result is then a tuple with the mean in the
            zeroth entry and the variance in the first entry. (default: False)

        """
Ultima's avatar
Ultima committed
180
181
        # Case 1: no operator given. Check function and domain for general
        # sanity
Ultima's avatar
Ultima committed
182
        if operator is None:
Ultima's avatar
Ultima committed
183
184
            # check whether the given function callable
            if function is None or not hasattr(function, "__call__"):
Ultima's avatar
Ultima committed
185
                raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
186
187
188
                  "ERROR: invalid input: No function given or not callable."))
            # check given domain
            if domain is None or not isinstance(domain, space):
Ultima's avatar
Ultima committed
189
                raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
190
                    "ERROR: invalid input: given domain is not a nifty space"))
191

Ultima's avatar
Ultima committed
192
193
        # Case 2: An operator is given. Take domain and function from that
        # if not given explicitly
Ultima's avatar
Ultima committed
194
        else:
195

Ultima's avatar
Ultima committed
196
197
198
            # Check 2.1 extract function
            # explicit function overrides operator function
            if function is None or not hasattr(function, "__call__"):
Ultima's avatar
Ultima committed
199
                try:
200
                    function = operator.times
Ultima's avatar
Ultima committed
201
202
                except(AttributeError):
                    raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
203
204
205
206
207
                        "ERROR: no explicit function given and given " +
                        "operator has no times method!"))

            # Check 2.2 check whether the given function is correctly bound to
            # the operator
Ultima's avatar
Ultima committed
208
209
            if operator != function.im_self:
                    raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
                        "ERROR: the given function is not a bound function " +
                        "of the operator!"))

            # Check 2.3 extract domain
            if domain is None or not isinstance(domain, space):
                if (function in [operator.inverse_times,
                                 operator.adjoint_times]):
                    try:
                        domain = operator.target
                    except(AttributeError):
                        raise ValueError(about._errors.cstring(
                            "ERROR: no explicit domain given and given " +
                            "operator has no target!"))
                else:
                    try:
                        domain = operator.domain
                    except(AttributeError):
                        raise ValueError(about._errors.cstring(
                            "ERROR: no explicit domain given and given " +
                            "operator has no domain!"))
230
231

        self.function = function
Ultima's avatar
Ultima committed
232
233
        self.domain = domain

Ultima's avatar
Ultima committed
234
        # Check the given codomain
Ultima's avatar
Ultima committed
235
236
        if codomain is None:
            codomain = self.domain.get_codomain()
Ultimanet's avatar
Ultimanet committed
237
        else:
Ultima's avatar
Ultima committed
238
            assert(self.domain.check_codomain(codomain))
239

Ultima's avatar
Ultima committed
240
241
        self.codomain = codomain

Ultima's avatar
Ultima committed
242
        if random not in ["pm1", "gau"]:
Ultima's avatar
Ultima committed
243
            raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
244
                "ERROR: unsupported random key '" + str(random) + "'."))
Ultima's avatar
Ultima committed
245
        self.random = random
246

Ultima's avatar
Ultima committed
247
        # Parse the remaining arguments
248
        self.nrun = int(nrun)
Ultima's avatar
Ultima committed
249
250
251
252
        self.varQ = bool(varQ)
        self.kwargs = kwargs

    def configure(self, **kwargs):
Ultimanet's avatar
Ultimanet committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
        """
            changes the attributes of the instance

            Parameters
            ----------
            random : string, *optional*
                the random number generator used to create the probes (default: "pm1")
            ncpu : int, *optional*
                the number of cpus to be used for parallel probing. (default: 2)
            nrun : int, *optional*
                the number of probes to be evaluated. If `nrun<ncpu`, it will be
                set to `ncpu`. (default: 8)
            nper : int, *optional*
                number of probes, that will be evaluated by one worker (default: 1)
            var : bool, *optional*
                whether the variance will be additionally returned (default: False)

        """
Ultima's avatar
Ultima committed
271
272
        if "random" in kwargs:
            if kwargs.get("random") not in ["pm1", "gau"]:
Ultima's avatar
Ultima committed
273
                raise ValueError(about._errors.cstring(
Ultima's avatar
Ultima committed
274
275
                    "ERROR: unsupported random key '" +
                    str(kwargs.get("random")) + "'."))
Ultimanet's avatar
Ultimanet committed
276
            else:
Ultima's avatar
Ultima committed
277
278
279
280
                self.random = kwargs.get("random")

        if "nrun" in kwargs:
            self.nrun = int(kwargs.get("nrun"))
Ultimanet's avatar
Ultimanet committed
281

Ultima's avatar
Ultima committed
282
283
        if "varQ" in kwargs:
            self.varQ = bool(kwargs.get("varQ"))
Ultimanet's avatar
Ultimanet committed
284

Ultima's avatar
Ultima committed
285
    def generate_probe(self):
Ultimanet's avatar
Ultimanet committed
286
287
288
289
290
291
292
293
294
295
        """
            Generates a single probe

            Returns
            -------
            probe : field
                a random field living in `domain` with mean 0 and variance 1 in
                each component

        """
Ultima's avatar
Ultima committed
296
        return field(self.domain,
Ultima's avatar
Ultima committed
297
298
                     codomain=self.codomain,
                     random=self.random)
Ultimanet's avatar
Ultimanet committed
299

Ultima's avatar
Ultima committed
300
    def evaluate_probe(self, probe, idnum=0):
Ultimanet's avatar
Ultimanet committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
        """
            Computes a single probing result given one probe

            Parameters
            ----------
            probe : field
                the field on which `function` will be applied
            idnum : int
                    the identification number of the probing

            Returns
            -------
            result : array-like
                the result of applying `function` to `probe`. The exact type
                depends on the function.

        """
Ultima's avatar
Ultima committed
318
319
        f = self.function(probe, **self.kwargs)
        return f
320

Ultima's avatar
Ultima committed
321
    def finalize(self, sum_of_probes, sum_of_squares, num):
Ultimanet's avatar
Ultimanet committed
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
        """
            Evaluates the probing results.

            Parameters
            ----------
            summa : numpy.array
                Sum of all probing results.
            num : int
                Number of successful probings (not returning ``None``).
            var : numpy.array
                Sum of all squared probing results

            Returns
            -------
            final : numpy.array
                The final probing result; 1/N Sum[ probing(probe_i) ]
            var : array-like
                The variance of the final probing result;
                1/(N(N-1)) Sum[( probing(probe_i) - final )^2];
                if the variance is returned, the return will be a tuple of
                (`final`,`var`).

        """
Ultima's avatar
Ultima committed
345
        # Check the success and efficiency of the probing
Ultima's avatar
Ultima committed
346
347
348
349
350
        if num < self.nrun:
            about.infos.cflush(
            " ( %u probe(s) failed, effectiveness == %.1f%% )\n"\
                %(self.nrun-num, 100*num/self.nrun))
            if num == 0:
Ultimanet's avatar
Ultimanet committed
351
352
353
354
355
                about.warnings.cprint("WARNING: probing failed.")
                return None
        else:
            about.infos.cflush("\n")

356

Ultima's avatar
Ultima committed
357
358
359
360
361
        if self.varQ == True:
            if num == 1:
                about.warnings.cprint(
                "WARNING: Only one probe available -> infinite variance.")
                return (sum_of_probes, None)
Ultimanet's avatar
Ultimanet committed
362
            else:
Ultima's avatar
Ultima committed
363
364
                var = 1/(num-1)*(sum_of_squares - 1/num*(sum_of_probes**2))
                return (sum_of_probes*(1./num), var)
Ultimanet's avatar
Ultimanet committed
365
        else:
Ultima's avatar
Ultima committed
366
            return sum_of_probes*(1./num)
Ultimanet's avatar
Ultimanet committed
367
368


Ultima's avatar
Ultima committed
369
370

    def print_progress(self, num): # > prints progress status by in upto 10 dots
Ultima's avatar
Ultima committed
371
        tenths = 1+(10*num//self.nrun)
Ultimanet's avatar
Ultimanet committed
372
        about.infos.cflush(("\b")*10+('.')*tenths+(' ')*(10-tenths))
Ultima's avatar
Ultima committed
373
    """
Ultima's avatar
Ultima committed
374
375
    def _single_probing(self,zipped): # > performs one probing operation
        # generate probe
Ultimanet's avatar
Ultimanet committed
376
377
        np.random.seed(zipped[0])
        probe = self.gen_probe()
Ultima's avatar
Ultima committed
378
        # do the actual probing
Ultimanet's avatar
Ultimanet committed
379
        return self.probing(zipped[1],probe)
Ultima's avatar
Ultima committed
380
    """
Ultimanet's avatar
Ultimanet committed
381

Ultima's avatar
Ultima committed
382
383
    def probe(self): # > performs the probing operations one after another
        # initialize the variables
Ultima's avatar
Ultima committed
384
385
386
        sum_of_probes = 0
        sum_of_squares = 0
        num = 0
387

Ultimanet's avatar
Ultimanet committed
388
        for ii in xrange(self.nrun):
389
            print ('running probe ', ii)
Ultima's avatar
Ultima committed
390
            temp_probe = self.generate_probe()
Ultima's avatar
Ultima committed
391
            temp_result = self.evaluate_probe(probe=temp_probe)
392

Ultima's avatar
Ultima committed
393
394
            if temp_result is not None:
                sum_of_probes += temp_result
395
396
                if self.varQ:
                    sum_of_squares += ((temp_result).conjugate())*temp_result
Ultima's avatar
Ultima committed
397
398
                num += 1
                self.print_progress(num)
Ultimanet's avatar
Ultimanet committed
399
        about.infos.cflush(" done.")
Ultima's avatar
Ultima committed
400
        # evaluate
Ultima's avatar
Ultima committed
401
        return self.finalize(sum_of_probes, sum_of_squares, num)
Ultimanet's avatar
Ultimanet committed
402

Ultima's avatar
Ultima committed
403
    def __call__(self, loop=False, **kwargs):
Ultimanet's avatar
Ultimanet committed
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
        """

            Starts the probing process.
            All keyword arguments that can be given to `configure` can also be
            given to `__call__` and have the same effect.

            Parameters
            ----------
            loop : bool, *optional*
                if `loop` is True, then multiprocessing will be disabled and
                all probes are evaluated by a single worker (default: False)

            Returns
            -------
            results : see **Returns** in `evaluate`

            other parameters
            ----------------
            kwargs : see **Parameters** in `configure`

        """
        self.configure(**kwargs)
Ultima's avatar
Ultima committed
426
        return self.probe()
Ultimanet's avatar
Ultimanet committed
427
428
429
430

    def __repr__(self):
        return "<nifty_core.probing>"

Ultima's avatar
Ultima committed
431

Ultima's avatar
Ultima committed
432
433
class _specialized_prober(object):
    def __init__(self, operator, domain=None, inverseQ=False, **kwargs):
Ultima's avatar
Ultima committed
434
        # remove a potentially supplied function keyword argument
Ultimanet's avatar
Ultimanet committed
435
        try:
Ultima's avatar
Ultima committed
436
437
438
            kwargs.pop('function')
        except(KeyError):
            pass
Ultimanet's avatar
Ultimanet committed
439
        else:
Ultima's avatar
Ultima committed
440
441
            about.warnings.cprint(
                "WARNING: Dropped the supplied function keyword-argument!")
442

Ultima's avatar
Ultima committed
443
        if domain is None and not inverseQ:
Ultima's avatar
Ultima committed
444
            kwargs.update({'domain': operator.domain})
Ultima's avatar
Ultima committed
445
        elif domain is None and inverseQ:
Ultima's avatar
Ultima committed
446
            kwargs.update({'domain': operator.target})
Ultimanet's avatar
Ultimanet committed
447
        else:
Ultima's avatar
Ultima committed
448
            kwargs.update({'domain': domain})
449
450
        self.operator = operator

Ultima's avatar
Ultima committed
451
        self.prober = prober(function=self._probing_function,
Ultima's avatar
Ultima committed
452
                             **kwargs)
453

Ultima's avatar
Ultima committed
454
455
    def _probing_function(self, probe):
        return None
456

Ultima's avatar
Ultima committed
457
458
    def __call__(self, *args, **kwargs):
        return self.prober(*args, **kwargs)
459

Ultima's avatar
Ultima committed
460
461
    def __getattr__(self, attr):
        return getattr(self.prober, attr)
462
463


Ultima's avatar
Ultima committed
464
class trace_prober(_specialized_prober):
465
    def __init__(self, operator, **kwargs):
Ultima's avatar
Ultima committed
466
        super(trace_prober, self).__init__(operator=operator,
467
                                           inverseQ=False,
Ultima's avatar
Ultima committed
468
                                           **kwargs)
Ultima's avatar
Ultima committed
469

Ultima's avatar
Ultima committed
470
    def _probing_function(self, probe):
471
        return direct_vdot(probe.conjugate(), self.operator.times(probe))
Ultima's avatar
Ultima committed
472

Ultima's avatar
Ultima committed
473

Ultima's avatar
Ultima committed
474
class inverse_trace_prober(_specialized_prober):
475
    def __init__(self, operator, **kwargs):
Ultima's avatar
Ultima committed
476
        super(inverse_trace_prober, self).__init__(operator=operator,
477
                                                   inverseQ=True,
Ultima's avatar
Ultima committed
478
                                                   **kwargs)
Ultima's avatar
Ultima committed
479

Ultima's avatar
Ultima committed
480
    def _probing_function(self, probe):
481
        return direct_vdot(probe.conjugate(),
Ultima's avatar
Ultima committed
482
483
                          self.operator.inverse_times(probe))

Ultima's avatar
Ultima committed
484

Ultima's avatar
Ultima committed
485
class diagonal_prober(_specialized_prober):
486
    def __init__(self, **kwargs):
Ultima's avatar
Ultima committed
487
        super(diagonal_prober, self).__init__(inverseQ=False,
Ultima's avatar
Ultima committed
488
                                              **kwargs)
Ultima's avatar
Ultima committed
489

Ultima's avatar
Ultima committed
490
491
    def _probing_function(self, probe):
        return probe.conjugate()*self.operator.times(probe)
492

Ultima's avatar
Ultima committed
493

Ultima's avatar
Ultima committed
494
class inverse_diagonal_prober(_specialized_prober):
495
    def __init__(self, operator, **kwargs):
Ultima's avatar
Ultima committed
496
        super(inverse_diagonal_prober, self).__init__(operator=operator,
497
                                                      inverseQ=True,
Ultima's avatar
Ultima committed
498
                                                      **kwargs)
Ultima's avatar
Ultima committed
499

Ultima's avatar
Ultima committed
500
501
    def _probing_function(self, probe):
        return probe.conjugate()*self.operator.inverse_times(probe)