diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e1a16d1df954983228b3b9778eb6fe9015b0a629..4a13e1831363718440c95b368c18ef6fc9e377b0 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -52,19 +52,23 @@ upload-job: - apt-get install -y -qq curl jq --fix-missing - apt-get install -y -qq python3-dev python3-pip # get the version that is being installed - - 'export version=$(grep -oP "(?<=__version__ = ).*" ./pymagglobal/__init__.py)' - - 'export version=${version:1:-1}' - - 'echo $version' + - export version=$(grep -oP "(?<=__version__ = ).*" ./pymagglobal/__init__.py) + - export version=${version:1:-1} + - echo $version # get the corresponding package id - - 'export ID=$(curl -s --request GET --header "PRIVATE-TOKEN: ${API_ACCES_TOKEN}" "https://git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages" | jq --arg VERSION "$version" --raw-output ".[] | if .version==$VERSION then .id else null end" | grep -v null)' - - 'echo $ID' + - | + export ID=$(curl -s --request GET --header "PRIVATE-TOKEN: ${API_ACCES_TOKEN}" "https://git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages" | jq --arg VERSION "$version" --raw-output '.[] | if .version==$VERSION then .id else null end' | grep -v null) + - echo $ID # remove the package, as we want a new build even if the version number did not change - - 'curl --max-time 900 --connect-timeout 120 --request DELETE --header "PRIVATE-TOKEN: ${API_ACCES_TOKEN}" "https://git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages/$ID"' + - | + curl --max-time 900 --connect-timeout 120 --request DELETE --header "PRIVATE-TOKEN: ${API_ACCES_TOKEN}" "https://git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages/$ID" # install twine - pip install twine # upload - TWINE_PASSWORD=${CI_JOB_TOKEN} TWINE_USERNAME=gitlab-ci-token python3 -m twine upload --repository-url https://git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages/pypi dist/* --verbose only: + refs: + - master changes: - pymagglobal/*.py - pymagglobal/dat/* @@ -81,6 +85,8 @@ test-conda: - pip install pymagglobal[tests] --extra-index-url https://public:5mz_iyigu-WE3HySBH1J@git.gfz-potsdam.de/api/v4/projects/${CI_PROJECT_ID}/packages/pypi/simple - python3 tests/run_tests.py only: + refs: + - master changes: - pymagglobal/*.py - pymagglobal/dat/* @@ -105,6 +111,14 @@ test-debian: - ./docs/pic_dip.png - ./docs/pic_cfe.png - ./docs/pic_cfs.png + only: + refs: + - master + changes: + - pymagglobal/*.py + - pymagglobal/dat/* + - setup.py + - tests/* pages: stage: deploy diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 index 0000000000000000000000000000000000000000..00dbef8449234572ede5444577f9e03a397da2d2 --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,4 @@ +2021-05-12 v0.1.1 + * the python backend can now handle models with uncertainty + * the field method now handles the field type (dif vs nez.) + * minor bug-fixes diff --git a/pymagglobal/_commands.py b/pymagglobal/_commands.py index 6fc584603446e0c66af5cde0b8390f26a9129d18..91a301e8e9ed12af051906b57600a2dff9923331 100644 --- a/pymagglobal/_commands.py +++ b/pymagglobal/_commands.py @@ -390,10 +390,7 @@ def maps(args): # BUGFIX: white edges in cartopy z_at[1] -= 180 # use the core function to get the field - field = core.field(z_at, args.model) - # convert the field if necessary - if args.type == 'dif': - field = np.array(utils.nez2dif(*field)) + field = core.field(z_at, args.model, field_type=args.type) # output formats for dif and nez components fmts = {'dif': ('%.1f', '%.1f', '%2.6f', '%2.6f', '%1.7e'), 'nez': ('%.1f', '%.1f', '%1.7e', '%1.7e', '%1.7e')} diff --git a/pymagglobal/core.py b/pymagglobal/core.py index 4f94de9a61f8910abdc5e41849e7a0f51d665137..1ebd76ef547700c55f86803b1bbb4558cb49a871 100644 --- a/pymagglobal/core.py +++ b/pymagglobal/core.py @@ -27,7 +27,7 @@ import numpy as np from scipy.interpolate import BSpline from pymagglobal.utils import nez2dif, REARTH, i2lm_l, i2lm_m, lmax2N, lm2i, \ - dsh_basis + dsh_basis, grad_d, grad_i, grad_f __all__ = ['built_in_models', 'Model', 'master_curve', 'coefficients', @@ -79,6 +79,9 @@ class Model(object): A large array containing coefficients at the (inner) knots. splines : scipy.interpolate.BSpline A `BSpline` instance representing the model. + cov_splines : scipy.interpolate.BSpline + A `BSPline` instance giving the covariances. Only available for certain + models. Default is `None`. ''' def __init__(self, name, fname=None): self.name = name @@ -107,6 +110,7 @@ class Model(object): self.splines = BSpline(self.knots, self.coeffs, 3) + self.cov_splines = None def valid_epoch(self, epoch): '''Check whether the given epoch is in the range of the model. @@ -124,12 +128,21 @@ class Model(object): epoch = np.asarray(epoch) out_of_range = epoch[(epoch < self.t_min) + (self.t_max < epoch)] if np.any(out_of_range): - warn(f'Epochs {out_of_range} are outside of the models range [{self.t_min},' - f' {self.t_max}]. The result will be extrapolated and not ' - f'robust.', + warn(f'Epochs {out_of_range} are outside of the models range ' + f'[{self.t_min}, {self.t_max}]. The result will be ' + f'extrapolated and not robust.', category=OutOfRangeWarning) return epoch + @property + def _has_covariances(self): + return (self.cov_splines is not None) + + @_has_covariances.setter + def _has_covariances(self, value): + raise RuntimeError("Private property '_has_covariances' can not be " + " set manually.") + def valid_degrees_orders(self, degrees, orders): '''Check whether the given degrees and orders are valid for the model and calculate the corresponding indices. @@ -163,7 +176,8 @@ class Model(object): return inds -def master_curve(times, loc, splines, ser_type='dif', check=True): +def master_curve(times, loc, splines, cov_splines=None, ser_type='dif', + check=True): '''Create master curves from a splines object. Parameters @@ -177,6 +191,10 @@ def master_curve(times, loc, splines, ser_type='dif', check=True): ser_type : {'dif', 'nez'} The type of the master curves. May be either 'dif' (default) for declination, inclination and intensity or 'nez' for north, east, down. + cov_splines : scipy.interpolate.BSpline, optional + A BSpline interpolating the covariance matrices of the coefficients. + If a model is passed to splines, this argument will be ignored and set, + if the model's covariances are available. check : bool, optional If a Model is given, check the input validity. @@ -184,18 +202,26 @@ def master_curve(times, loc, splines, ser_type='dif', check=True): ------- ndarray The first component master curve. Either declination or north, - depending on the ser_type kwarg. + depending on the ser_type kwarg. If cov_splines is given, a tuple + containing the component and the standard deviation is returned. ndarray The second component master curve. Either inclination or east, - depending on the ser_type kwarg. + depending on the ser_type kwarg. If cov_splines is given, a tuple + containing the component and the standard deviation is returned. ndarray The third component master curve. Either intensity or down, - depending on the ser_type kwarg. + depending on the ser_type kwarg. If cov_splines is given, a tuple + containing the component and the standard deviation is returned. ''' + if ser_type not in ['nez', 'dif']: + raise ValueError(f"{ser_type} is not a valid series type. 'dif' and " + f"'nez' are supported.") if isinstance(splines, Model): if check: for it in np.atleast_1d(times): splines.valid_epoch(it) + if splines._has_covariances: + cov_splines = splines.cov_splines splines = splines.splines try: n_plt = len(times) @@ -208,27 +234,31 @@ def master_curve(times, loc, splines, ser_type='dif', check=True): z_mst[2] = REARTH z_mst[3] = times # the master curves are just field values at a constant location - mst_cv = field(z_mst, splines) + mst_cv = field(z_mst, splines, cov_splines=cov_splines, + field_type=ser_type) - if ser_type == 'dif': - return nez2dif(*mst_cv) - elif ser_type == 'nez': + if cov_splines is None: return mst_cv[0], mst_cv[1], mst_cv[2] - else: - raise ValueError(f"{ser_type} is not a valid series type. 'dif' and " - f"'nez' are supported.") + return (mst_cv[0][0], mst_cv[1][0]), (mst_cv[0][1], mst_cv[1][1]), \ + (mst_cv[0][2], mst_cv[1][2]) -def coefficients(epoch, splines, check=True): + +def coefficients(epoch, splines, cov_splines=None, check=True): '''Evaluate splines at an epoch and return the coefficients, together with the appropriate degrees and orders. Parameters ---------- - epoch : float - The epoch at which to return the coefficients, in years. + epoch : float or array + The epoch(s) at which to return the coefficients, in years. splines : scipy.interpolate.BSpline or Model An instance of Model or splines specifying the model. + cov_splines : scipy.interpolate.BSpline, optional + A BSpline interpolating the covariance matrices of the coefficients. + If a model is passed to splines, this argument will be ignored and set, + if the model's covariances are available. If cov_splines is given, + the standard deviations are returned as a fourth return value. check : bool, optional If a Model is given, check the input validity. @@ -246,16 +276,27 @@ def coefficients(epoch, splines, check=True): if check: for it in np.atleast_1d(epoch): splines.valid_epoch(epoch) + if splines._has_covariances: + cov_splines = splines.cov_splines splines = splines.splines + coeffs = splines(epoch) # set up degrees and orders for output purposes ls = [i2lm_l(it) for it in range(splines.c.shape[1])] ms = [i2lm_m(it) for it in range(splines.c.shape[1])] + if cov_splines is None: + return ls, ms, coeffs - return ls, ms, coeffs + covs = cov_splines(epoch) + try: + covs = np.diagonal(covs, axis1=1, axis2=2) + except np.AxisError: + covs = np.diagonal(covs) + return ls, ms, coeffs, np.sqrt(covs) -def dipole_series(times, splines, check=True): + +def dipole_series(times, splines, cov_splines=None, check=True): '''Create a dipole-moment time series from a splines object. Parameters @@ -264,28 +305,54 @@ def dipole_series(times, splines, check=True): The times for which to create the master curve. splines : scipy.interpolate.BSpline or Model An instance of Model or splines specifying the model. + cov_splines : scipy.interpolate.BSpline, optional + A BSpline interpolating the covariance matrices of the coefficients. + If a model is passed to splines, this argument will be ignored and set, + if the model's covariances are available. If cov_splines is given, + the 16- and 84-percentiles are returned along with the series. check : bool, optional If a Model is given, check the input validity. Returns ------- ndarray - An array containing the dipole-moment time series. + An array containing the dipole-moment time series. If cov_splines is + given, the 16- and 84-percentiles are returned along with the series. ''' if isinstance(splines, Model): if check: for it in np.atleast_1d(times): splines.valid_epoch(it) + if splines._has_covariances: + cov_splines = splines.cov_splines splines = splines.splines - # evaluate the splines - coeffs = splines(times).T - # calculate the dipole moment - dip_ser = np.sqrt(coeffs[0]**2 + coeffs[1]**2 - + coeffs[2]**2) - # get the units to 1e22 [m^2 A] - dip_ser *= 0.63712**3/1000 - return dip_ser + # conversion factor from nT to 1e22 [m^2 A] + conversion = 0.63712**3/1000 + if cov_splines is None: + # evaluate the splines + coeffs = splines(times).T + # calculate the dipole moment + dip_ser = np.sqrt(coeffs[0]**2 + coeffs[1]**2 + + coeffs[2]**2) + dip_ser *= conversion + return dip_ser + + coeffs = splines(times)[:, :3] + covs = cov_splines(times)[:, :3, :3] + mus = np.empty(len(times)) + errs_16 = np.empty(len(times)) + errs_84 = np.empty(len(times)) + for it, (mu, cov) in enumerate(zip(coeffs, covs)): + ens_dip = np.random.multivariate_normal(mean=mu, + cov=cov, + size=10000).T + ens_dip = ens_dip.reshape(-1, 3, 10000) + ens_m = np.sqrt(np.square(ens_dip).sum(axis=1)) + ens_m *= conversion + mus[it] = ens_m.mean(axis=1) + errs_16[it], errs_84[it] = np.percentile(ens_m, (16, 84), axis=1) + return mus, errs_16, errs_84 def file2splines(fname): @@ -312,9 +379,8 @@ def file2splines(fname): return c_splines -def field(z_at, splines, check=True): - '''Evaluate coefficient splines at locations z_at and return the field in - north, east, down components. +def field(z_at, splines, cov_splines=None, field_type='nez', check=True): + '''Evaluate coefficient splines at locations z_at and return the field. Parameters ---------- @@ -331,19 +397,30 @@ def field(z_at, splines, check=True): points. splines : scipy.interpolate.BSpline or Model An instance of Model or splines specifying the model. + cov_splines : scipy.interpolate.BSpline, optional + A BSpline interpolating the covariance matrices of the coefficients. + If a model is passed to splines, this argument will be ignored and set, + if the model's covariances are available. If cov_splines is given, + the standard deviations are returned as a second return value. + field_type : {'dif', 'nez'} + The type of output. May be either 'nez' (default) for north, east and + down component or 'dif' for declination, inclination and intensity. check : bool, optional If a Model is given, check the input validity. Returns ------- array of shape (3, n) - An array containing the north, east and down component at the input - locations. + An array containing the field components at the input locations. If + cov_splines is given, the standard deviations are returned as a + second return value of identical shape. ''' if isinstance(splines, Model): if check: for it in np.atleast_1d(z_at[3]): splines.valid_epoch(it) + if splines._has_covariances: + cov_splines = splines.cov_splines splines = splines.splines # count the number of inputs n_pts = z_at.shape[1] @@ -361,6 +438,45 @@ def field(z_at, splines, check=True): # calculate the field as a dot product between coefficients and basis field = np.einsum('ij, ijk -> jk', coeffs, - spatial) - - return field.T + spatial).T + + if cov_splines is None: + if field_type == 'dif': + field = np.array(nez2dif(*field)) + return field + + covs = cov_splines(z_at[3]).T + + if field_type == 'dif': + field_cov = np.einsum('ilk, ijl, jlm -> lkm', + spatial, + covs, + spatial) + + grad_D = grad_d(*field) + cov_D = np.einsum('ij, ijk, ik -> i', + grad_D, + field_cov, + grad_D) + + grad_I = grad_i(*field) + cov_I = np.einsum('ij, ijk, ik -> i', + grad_I, + field_cov, + grad_I) + + grad_F = grad_f(*field) + cov_F = np.einsum('ij, ijk, ik -> i', + grad_F, + field_cov, + grad_F) + + return np.array(nez2dif(*field)), \ + np.sqrt(np.array([cov_D, cov_I, cov_F])) + + field_cov = np.einsum('ilk, ijl, jlk -> lk', + spatial, + covs, + spatial) + + return field, np.sqrt(field_cov).T diff --git a/pymagglobal/utils.py b/pymagglobal/utils.py index 89d9a3caf45b40829d48269090933d3434f8dc74..fab697974eb4251c853c95a5442fec2423318095 100644 --- a/pymagglobal/utils.py +++ b/pymagglobal/utils.py @@ -39,6 +39,32 @@ def _nicelabel(label): return fr'${parts[0]}$ {parts[1]}' +def grad_d(n, e, z): + '''Calculate the gradient of D. + ''' + res = np.column_stack((-e, n, np.zeros_like(n))) + res /= (n**2 + e**2)[:, np.newaxis] + return np.rad2deg(res) + + +def grad_i(n, e, z): + '''Calculate the gradient of I. + ''' + res = -np.column_stack((n, e, z)) \ + * z[:, np.newaxis] + res /= (n**2 + e**2 + z**2)[:, np.newaxis] + res[:, 2] += 1. + res /= np.sqrt(n**2 + e**2)[:, np.newaxis] + return np.rad2deg(res) + + +def grad_f(n, e, z): + '''Calculate the gradient of F ''' + res = np.column_stack((n, e, z)) + res /= np.sqrt(n**2 + e**2 + z**2)[:, np.newaxis] + return res + + def nez2dif(n, e, z): '''Transform the magnetic field components north, east and vertically downward to declination, inclination and intensity.