Commit fb536d3b authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

License note update.

parent 00b28c80
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This file is part of CORBASS.
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +19,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
from . import utils
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This module is part of the CORBASS algorithm. It provides routines for the
third and final step, the evaluation of the output from the integration
step. These routines serve to calculate probability density functions (or
proxys of these) for quantities of interest, such as the power spectrum,
the dipole intensity, the dipole location etc.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +25,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
import sys
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This module is part of the CORBASS algorithm. It provides a routine for the
first step, the exploration of the parameter space.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +22,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
import sys
......@@ -34,6 +37,7 @@ import numpy as np
import progressbar
from corbass.inversion import Inversion
from corbass.utils import bounds2grid
class ExplorationResult:
......@@ -104,29 +108,16 @@ def explore(name, data, N, bounds, r_ref=3200., prefix=''):
"""
# set up exploration grid
lam_arr, dlam = np.linspace(bounds[0, 0],
bounds[0, 1],
N,
retstep=True)
eps_arr, deps = np.linspace(bounds[1, 0],
bounds[1, 1],
N,
retstep=True)
rho_arr, drho = np.linspace(bounds[2, 0],
bounds[2, 1],
N,
retstep=True)
lams, epss, rhos = bounds2grid(bounds, N)
# set up inversion class for the data
# @Max: What happens if a bin is empty?
inv = Inversion(data, drop_duplicates=False)
# set up iterator
it = np.nditer([lam_arr[:, None, None],
eps_arr[None, :, None],
rho_arr[None, None, :],
it = np.nditer([lams[0][:, None, None],
epss[0][None, :, None],
rhos[0][None, None, :],
None])
# set up progressbar
......@@ -149,7 +140,7 @@ def explore(name, data, N, bounds, r_ref=3200., prefix=''):
# log-posterior to posterior
posterior = np.exp(posterior)
# normalize
posterior /= dlam*deps*drho * posterior.sum()
posterior /= lams[1]*epss[1]*rhos[1] * posterior.sum()
# save results
rslt = ExplorationResult(name)
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This module is part of the CORBASS algorithm. It provides routines for the
second step, the integration of a section of the parameter space. These
consist of a routine for the actual integration step, as well as an
assisting routine for finding a region of high probability mass, given
output from the previous exploration step.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +25,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
import sys
......@@ -33,6 +39,7 @@ import numpy as np
import progressbar
from corbass.inversion import Inversion
from corbass.utils import bounds2grid
class IntegrationResult:
......@@ -161,43 +168,6 @@ class IntegrationResult:
print('No results available, nothing to write!')
def bounds2grid(expl_bounds, N):
""" Set up 3 arrays of length N from the exploration boundaries.
Parameters
----------
expl_bounds : 3x2 array-like
The boundaries from the exploration, as min-max pairs for the
non-dipole variance, the error scaling and the residual
N : int
The number of points per dimension
Returns
-------
lams : tuple
A tuple of the non-dipole variance array and the gridstep
epss : tuple
A tuple of the error scaling array and the gridstep
rhoss : tuple
A tuple of the residual array and the gridstep
"""
lams = np.linspace(expl_bounds[0, 0],
expl_bounds[0, 1],
N,
retstep=True)
epss = np.linspace(expl_bounds[1, 0],
expl_bounds[1, 1],
N,
retstep=True)
rhos = np.linspace(expl_bounds[2, 0],
expl_bounds[2, 1],
N,
retstep=True)
return lams, epss, rhos
def get_integ_bounds(expl_bounds, posterior, moments_out=False):
""" From the exploration grid and the posterior, calculate the first two
moments for each marginal distribution to set up bounds for the
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This module is part of the CORBASS algorithm. It can be seen as the core
module, where the actual inversion code is defined.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,9 +22,8 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
if __name__ == '__main__':
# make relative imports work
import sys
......@@ -35,8 +37,7 @@ import pyfield
from corbass.utils import scaling, read_data
r_Ear = pyfield.REARTH
r_CMB = 3480.
RCMB = 3480.
def grad_d(b):
......@@ -183,7 +184,7 @@ class Inversion:
A residual error.
"""
# Kernels from the pyfield library
self.sp_Dip = pyfield.StatDipIntSpline(r_Ear)
self.sp_Dip = pyfield.StatDipIntSpline(pyfield.REARTH)
self.sp_WoDip = pyfield.StatWODipIntSpline(r_Ref)
# Allocate and retrieve basis functions (-grad Phi)
self.dspharm_obs = np.empty((3, 3*self.n_obs), order='F')
......@@ -480,7 +481,7 @@ class Inversion:
the corresponding pointwise variance
"""
# Setup design points
x_desi, n_desi = self.desi_points(n_desi, rad=r_CMB)
x_desi, n_desi = self.desi_points(n_desi, rad=RCMB)
# Allocate and retrieve basis functions at design points
dspharm_desi = np.empty((3, 3*n_desi), order='F')
......@@ -687,7 +688,7 @@ class Inversion:
dummy))
# Dipole part from the non-informative prior
sss = scaling(r_Ear, r_Ref, 1)[0]
sss = scaling(pyfield.REARTH, r_Ref, 1)[0]
# Scale to R_ref
mu_g[0:3] = np.multiply(self.g_bar, sss)
cov_g[0:3, 0:3] = np.multiply(self.var_g_bar, sss**2)
......@@ -1058,7 +1059,7 @@ class Inversion:
self.idx_NEZ_I,
self.idx_NEZ_F))
def desi_points(self, n, rad=r_Ear, ret_n=True):
def desi_points(self, n, rad=pyfield.REARTH, ret_n=True):
""" Return equidistant design points on the sphere, as expected by
the pyfield library.
......@@ -1126,7 +1127,7 @@ if __name__ == '__main__':
alpha_WoDip = 16000
sigma_E = 1.34
sigma_R = 5000
x_zf, zs, fs = test.zf_CMB_inversion(r_ref, alpha_WoDip, sigma_E, sigma_R,
1000)
lat_zf = 90-x_zf[0]
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This module is part of the CORBASS algorithm. It provides utility routines
that are used throughut the algorithm, e.g. a routine to load a parameter
file or scaling routines for coefficicients and covariances.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +23,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
import sys
......@@ -98,6 +102,43 @@ def read_data(fname):
'D', 'dD', 'I', 'dI', 'F', 'dF']]
def bounds2grid(bounds, N):
""" Set up 3 arrays of length N from boundaries.
Parameters
----------
bounds : 3x2 array-like
The boundaries of a region in the parameter space, as min-max pairs
for the non-dipole variance, the error scaling and the residual
N : int
The number of points per dimension
Returns
-------
lams : tuple
A tuple of the non-dipole variance array and the gridstep
epss : tuple
A tuple of the error scaling array and the gridstep
rhoss : tuple
A tuple of the residual array and the gridstep
"""
lams = np.linspace(bounds[0, 0],
bounds[0, 1],
N,
retstep=True)
epss = np.linspace(bounds[1, 0],
bounds[1, 1],
N,
retstep=True)
rhos = np.linspace(bounds[2, 0],
bounds[2, 1],
N,
retstep=True)
return lams, epss, rhos
def newer(file_a, file_b):
""" Use the modification date to check if file_a is newer than file_b """
diff = Path(file_a).stat().st_mtime - Path(file_b).stat().st_mtime
......@@ -141,7 +182,6 @@ def load(argv):
except AttributeError:
final = data['t'].max()
# explicitly calculate bin edges
# +1 to include the final year
bins = range(initial, final + 1, pars.binlength)
......
__doc__ = """ This is an example parameter file for CORBASS. """
""" This is an example parameter file for CORBASS. """
from pathlib import Path
path = Path(__file__).parent # relative to the location of this very par file
......
"""
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
This script runs the CORBASS algorithm for a parameter file that is
passed via the command line.
This file is part of CORBASS.
Copyright (C) 2019 Maximilian Schanner and Stefan Mauerberger
Cite as:
Schanner, M. and Mauerberger, S. (2019)
CORBASS: CORrelation Based Archeomagnetic SnapShot model. V. 1.0.
GFZ Data Services. http://doi.org/XXX TODO XXX
This file is part of CORBASS.
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
......@@ -19,7 +22,7 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
along with CORBASS. If not, see <https://www.gnu.org/licenses/>.
"""
import sys
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment