exploration.py 4.81 KB
Newer Older
Maximilian Schanner's avatar
Maximilian Schanner committed
1
"""
Maximilian Schanner's avatar
Maximilian Schanner committed
2
3
    This module is part of the CORBASS algorithm. It provides a routine for the
    first step, the exploration of the parameter space.
Maximilian Schanner's avatar
Maximilian Schanner committed
4

Maximilian Schanner's avatar
Maximilian Schanner committed
5
6
    Copyright (C) 2019 Helmholtz Centre Potsdam GFZ,
        German Research Centre for Geosciences, Potsdam, Germany
Maximilian Schanner's avatar
Maximilian Schanner committed
7
8

    Cite as:
Maximilian Schanner's avatar
Maximilian Schanner committed
9
    Schanner, Maximilian Arthus and Mauerberger, Stefan (2019)
Maximilian Schanner's avatar
Maximilian Schanner committed
10
    CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
Maximilian Schanner's avatar
Maximilian Schanner committed
11
    GFZ Data Services. http://doi.org/10.5880/GFZ.2.3.2019.008
Maximilian Schanner's avatar
Maximilian Schanner committed
12

Maximilian Schanner's avatar
Maximilian Schanner committed
13
14
    This file is part of CORBASS.

Maximilian Schanner's avatar
Maximilian Schanner committed
15
16
17
18
19
20
21
22
23
24
25
    CORBASS 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.

    CORBASS 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
Maximilian Schanner's avatar
Maximilian Schanner committed
26
    along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
Maximilian Schanner's avatar
Maximilian Schanner committed
27
28
"""

Maximilian Schanner's avatar
Maximilian Schanner committed
29
30
31
32
import sys

if __name__ == '__main__':
    # make relative imports work
Maximilian Schanner's avatar
Maximilian Schanner committed
33
34
    from pathlib import Path
    sys.path.append(str(Path(__file__).absolute().parents[1]))
Maximilian Schanner's avatar
Maximilian Schanner committed
35
36
37
38
39

import numpy as np

import progressbar

Maximilian Schanner's avatar
Maximilian Schanner committed
40
from corbass.inversion import Inversion
Maximilian Schanner's avatar
Maximilian Schanner committed
41
from corbass.utils import bounds2grid
Maximilian Schanner's avatar
Maximilian Schanner committed
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


class ExplorationResult:
    """ Class for representing the results of an exploration run. During the
    run the Attributes will be set.

    Parameters
    ----------
    name : str
        A name identifying the result. This name will also specify the output
        name.

    Attributes
    ----------
    suffix : str
        The suffix appended to the output-filename.

    Methods
    -------
    write()
        Write the results to a file. The following attributes should be set by
        the routine returning an IntegrationResult:
            posterior : ndarray
                The posterior on a grid.
            scale : float
                A scaling for the posterior, which may be necessary if the
                values are small, so that np.exp(posterior) evaluates to 0. A
                good choice is the maximum of the posterior resulting from the
                grid exploration.
    """
    suffix = '_expl.npz'

    def __init__(self, name):
        self.name = name

    def write(self):
        try:
            np.savez(f'{self.name}{self.suffix}',
                     posterior=self.posterior,
                     scale=self.scale)
        except AttributeError:
            print('No results available, nothing to write!')


86
def explore(name, data, N, bounds, r_ref=3200.):
Maximilian Schanner's avatar
Maximilian Schanner committed
87
88
89
90
91
    """ Explore the parameterspace given by bounds on an NxNxN grid, for the
    given data and save the output

    Parameters
    ----------
92
93
    name : str
        A name for the bin
Maximilian Schanner's avatar
Maximilian Schanner committed
94
95
96
97
98
99
100
    data : Pandas.DataFrame
        The data to be explored
    N : int
        The number of points per dimension
    bounds : 3x2 array-like
        The boundaries for exploration, as min-max pairs for the non-dipole
        variance, the error scaling and the residual
101
    r_ref : float (optional, default is 3200.0)
Maximilian Schanner's avatar
Maximilian Schanner committed
102
103
104
105
106
107
108
109
110
111
        The reference radius

    Returns
    -------
    ExplorationResult
        The result from the exploration run as a specific object, containing
        the results as attributes and providing a write method.
    """

    # set up exploration grid
Maximilian Schanner's avatar
Maximilian Schanner committed
112
    lams, epss, rhos = bounds2grid(bounds, N)
Maximilian Schanner's avatar
Maximilian Schanner committed
113
114
115
116
117
118

    # set up inversion class for the data
    # @Max: What happens if a bin is empty?
    inv = Inversion(data, drop_duplicates=False)

    # set up iterator
Maximilian Schanner's avatar
Maximilian Schanner committed
119
120
121
    it = np.nditer([lams[0][:, None, None],
                    epss[0][None, :, None],
                    rhos[0][None, None, :],
Maximilian Schanner's avatar
Maximilian Schanner committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
                    None])

    # set up progressbar
    pbar = progressbar.ProgressBar(max_value=N**3,
                                   redirect_stdout=True)
    pbar.start()
    # iterate
    for lam, eps, rho, posterior in it:
        posterior[...] = inv.log_pst(r_ref, lam, eps, rho)

        pbar.update(pbar.value + 1)
    # close progressbar
    pbar.finish()

    # get posterior
    posterior = it.operands[-1]
    # rescale
    scale = posterior.max()
    posterior -= scale
    # log-posterior to posterior
    posterior = np.exp(posterior)
    # normalize
Maximilian Schanner's avatar
Maximilian Schanner committed
144
    posterior /= lams[1]*epss[1]*rhos[1] * posterior.sum()
Maximilian Schanner's avatar
Maximilian Schanner committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

    # save results
    rslt = ExplorationResult(name)
    rslt.posterior = posterior
    rslt.scale = scale

    return rslt


if __name__ == '__main__':
    from utils import load

    pars = load(sys.argv)

    for bin_name, data in pars.data:
        print(f"Exploring {bin_name}, this may take a while...")
        rslt = explore(pars.bin_fname(bin_name), data, pars.n_ex, pars.bounds,
                       r_ref=pars.r_ref)
        rslt.write()