476491407

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/GeoStat-Framework/GSTools/commit/<a class=hub.com/GeoStat-Framework/GSTools/commit/bb07527fd578bf23f96b516bd02bfece874db403">bb07527fd<a href="https://github.com/GeoStat-Framework/GSTools/commit/bb07527fd578bf23f96b516bd02bfece874db403">">Merge </a><a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/<a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/120756daf1243ac1a8ee313b55ad06af54f43fe9">120756daf</a>">120756daf</a><a href="https://github.com/GeoStat-Framework/GSTools/commit/bb07527fd578bf23f96b516bd02bfece874db403"> into e28b00c95">e28b00c95</a>

Pull Request Pull Request #113: Add support for geographical coordinates

Run Details

87 of 87 new or added lines in 8 files covered. (100.0%)

2101 of 2319 relevant lines covered (90.6%)

13.59 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.61

/gstools/covmodel/fit.py

# -*- coding: utf-8 -*-
"""
GStools subpackage providing tools for the covariance-model.

.. currentmodule:: gstools.covmodel.fit

The following classes and functions are provided

.. autosummary::
   fit_variogram
"""

# pylint: disable=C0103
import numpy as np
from scipy.optimize import curve_fit
from gstools.covmodel.tools import check_arg_in_bounds, default_arg_from_bounds
from gstools.tools.geometric import set_anis


__all__ = ["fit_variogram"]


DEFAULT_PARA = ["var", "len_scale", "nugget"]


def fit_variogram(
    model,
    x_data,
    y_data,
    anis=True,
    sill=None,
    init_guess="default",
    weights=None,
    method="trf",
    loss="soft_l1",
    max_eval=None,
    return_r2=False,
    curve_fit_kwargs=None,
    **para_select
):
    """
    Fitting a variogram-model to an empirical variogram.

    Parameters
    ----------
    model : :any:`CovModel`
        Covariance Model to fit.
    x_data : :class:`numpy.ndarray`
        The bin-centers of the empirical variogram.
    y_data : :class:`numpy.ndarray`
        The messured variogram
        If multiple are given, they are interpreted as the directional
        variograms along the main axis of the associated rotated
        coordinate system.
        Anisotropy ratios will be estimated in that case.
    anis : :class:`bool`, optional
        In case of a directional variogram, you can control anisotropy
        by this argument. Deselect the parameter from fitting, by setting
        it "False".
        You could also pass a fixed value to be set in the model.
        Then the anisotropy ratios won't be altered during fitting.
        Default: True
    sill : :class:`float` or :class:`bool` or :any:`None`, optional
        Here you can provide a fixed sill for the variogram.
        It needs to be in a fitting range for the var and nugget bounds.
        If variance or nugget are not selected for estimation,
        the nugget will be recalculated to fulfill:

            * sill = var + nugget
            * if the variance is bigger than the sill,
              nugget will bet set to its lower bound
              and the variance will be set to the fitting partial sill.

        If variance is deselected, it needs to be less than the sill,
        otherwise a ValueError comes up. Same for nugget.
        If sill=False, it will be deslected from estimation
        and set to the current sill of the model.
        Then, the procedure above is applied.
        Default: None
    init_guess : :class:`str`, optional
        Initial guess for the estimation. Either:

            * "default": using the default values of the covariance model
            * "current": using the current values of the covariance model

        Default: "default"
    weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`optional
        Weights applied to each point in the estimation. Either:

            * 'inv': inverse distance ``1 / (x_data + 1)``
            * list: weights given per bin
            * callable: function applied to x_data

        If callable, it must take a 1-d ndarray. Then ``weights = f(x_data)``.
        Default: None
    method : {'trf', 'dogbox'}, optional
        Algorithm to perform minimization.

            * 'trf' : Trust Region Reflective algorithm, particularly suitable
              for large sparse problems with bounds. Generally robust method.
            * 'dogbox' : dogleg algorithm with rectangular trust regions,
              typical use case is small problems with bounds. Not recommended
              for problems with rank-deficient Jacobian.

        Default: 'trf'
    loss : :class:`str` or :class:`callable`, optional
        Determines the loss function in scipys curve_fit.
        The following keyword values are allowed:

            * 'linear' (default) : ``rho(z) = z``. Gives a standard
              least-squares problem.
            * 'soft_l1' : ``rho(z) = 2 * ((1 + z)**0.5 - 1)``. The smooth
              approximation of l1 (absolute value) loss. Usually a good
              choice for robust least squares.
            * 'huber' : ``rho(z) = z if z <= 1 else 2*z**0.5 - 1``. Works
              similarly to 'soft_l1'.
            * 'cauchy' : ``rho(z) = ln(1 + z)``. Severely weakens outliers
              influence, but may cause difficulties in optimization process.
            * 'arctan' : ``rho(z) = arctan(z)``. Limits a maximum loss on
              a single residual, has properties similar to 'cauchy'.

        If callable, it must take a 1-d ndarray ``z=f**2`` and return an
        array_like with shape (3, m) where row 0 contains function values,
        row 1 contains first derivatives and row 2 contains second
        derivatives. Default: 'soft_l1'
    max_eval : :class:`int` or :any:`None`, optional
        Maximum number of function evaluations before the termination.
        If None (default), the value is chosen automatically: 100 * n.
    return_r2 : :class:`bool`, optional
        Whether to return the r2 score of the estimation.
        Default: False
    curve_fit_kwargs : :class:`dict`, optional
        Other keyword arguments passed to scipys curve_fit. Default: None
    **para_select
        You can deselect parameters from fitting, by setting
        them "False" using their names as keywords.
        You could also pass fixed values for each parameter.
        Then these values will be applied and the involved parameters wont
        be fitted.
        By default, all parameters are fitted.

    Returns
    -------
    fit_para : :class:`dict`
        Dictonary with the fitted parameter values
    pcov : :class:`numpy.ndarray`
        The estimated covariance of `popt` from
        :any:`scipy.optimize.curve_fit`.
        To compute one standard deviation errors
        on the parameters use ``perr = np.sqrt(np.diag(pcov))``.
    r2_score : :class:`float`, optional
        r2 score of the curve fitting results. Only if return_r2 is True.

    Notes
    -----
    You can set the bounds for each parameter by accessing
    :any:`CovModel.set_arg_bounds`.

    The fitted parameters will be instantly set in the model.
    """
    # preprocess selected parameters
    para, sill, constrain_sill, anis = _pre_para(
        model, para_select, sill, anis
    )
    # check curve_fit kwargs
    curve_fit_kwargs = {} if curve_fit_kwargs is None else curve_fit_kwargs
    # check method
    if method not in ["trf", "dogbox"]:
        raise ValueError("fit: method needs to be either 'trf' or 'dogbox'")
    # prepare variogram data
    # => concatenate directional variograms to have a 1D array for x and y
    x_data, y_data, is_dir_vario = _check_vario(model, x_data, y_data)
    # only fit anisotropy if a directional variogram was given
    anis &= is_dir_vario
    # set weights
    _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario)
    # set the lower/upper boundaries for the variogram-parameters
    bounds, init_guess_list = _init_curve_fit_para(
        model, para, init_guess, constrain_sill, sill, anis
    )
    # create the fitting curve
    curve_fit_kwargs["f"] = _get_curve(
        model, para, constrain_sill, sill, anis, is_dir_vario
    )
    # set the remaining kwargs for curve_fit
    curve_fit_kwargs["bounds"] = bounds
    curve_fit_kwargs["p0"] = init_guess_list
    curve_fit_kwargs["xdata"] = x_data
    curve_fit_kwargs["ydata"] = y_data
    curve_fit_kwargs["loss"] = loss
    curve_fit_kwargs["max_nfev"] = max_eval
    curve_fit_kwargs["method"] = method
    # fit the variogram
    popt, pcov = curve_fit(**curve_fit_kwargs)
    # convert the results
    fit_para = _post_fitting(model, para, popt, anis, is_dir_vario)
    # calculate the r2 score if wanted
    if return_r2:
        return fit_para, pcov, _r2_score(model, x_data, y_data, is_dir_vario)
    return fit_para, pcov


def _pre_para(model, para_select, sill, anis):
    """Preprocess selected parameters."""
    var_last = False
    for par in para_select:
        if par not in model.arg_bounds:
            raise ValueError(
                "fit: unknow parameter in selection: {}".format(par)
            )
        if not isinstance(para_select[par], bool):
            if par == "var":
                var_last = True
                var_tmp = float(para_select[par])
            else:
                setattr(model, par, float(para_select[par]))
            para_select[par] = False
    # set variance last due to possible recalculations
    if var_last:
        model.var = var_tmp
    # remove those that were set to True
    para_select = {k: v for k, v in para_select.items() if not v}
    # handling the sill
    sill = None if (isinstance(sill, bool) and sill) else sill
    if sill is not None:
        sill = model.sill if isinstance(sill, bool) else float(sill)
        constrain_sill = True
        sill_low = model.arg_bounds["var"][0] + model.arg_bounds["nugget"][0]
        sill_up = model.arg_bounds["var"][1] + model.arg_bounds["nugget"][1]
        if not sill_low <= sill <= sill_up:
            raise ValueError("fit: sill out of bounds.")
        if "var" in para_select and "nugget" in para_select:
            if model.var > sill:
                model.nugget = model.arg_bounds["nugget"][0]
                model.var = sill - model.nugget
            else:
                model.nugget = sill - model.var
        elif "var" in para_select:
            if model.var > sill:
                raise ValueError(
                    "fit: if sill is fixed and variance deselected, "
                    + "the set variance should be less than the given sill."
                )
            para_select["nugget"] = False
            model.nugget = sill - model.var
        elif "nugget" in para_select:
            if model.nugget > sill:
                raise ValueError(
                    "fit: if sill is fixed and nugget deselected, "
                    + "the set nugget should be less than the given sill."
                )
            para_select["var"] = False
            model.var = sill - model.nugget
        else:
            # deselect the nugget, to recalculate it accordingly
            # nugget = sill - var
            para_select["nugget"] = False
    else:
        constrain_sill = False
    # select all parameters to be fitted
    para = {par: True for par in DEFAULT_PARA}
    para.update({opt: True for opt in model.opt_arg})
    # now deselect unwanted parameters
    para.update(para_select)
    # check if anisotropy should be fitted or set
    if not isinstance(anis, bool):
        model.anis = anis
        anis = False
    return para, sill, constrain_sill, anis


def _check_vario(model, x_data, y_data):
    # prepare variogram data
    x_data = np.array(x_data).reshape(-1)
    y_data = np.array(y_data).reshape(-1)
    # if multiple variograms are given, they will be interpreted
    # as directional variograms along the main rotated axes of the model
    is_dir_vario = False
    if model.dim > 1 and x_data.size * model.dim == y_data.size:
        is_dir_vario = True
        # concatenate multiple variograms
        x_data = np.tile(x_data, model.dim)
    elif x_data.size != y_data.size:
        raise ValueError(
            "CovModel.fit_variogram: Wrong number of empirical variograms! "
            + "Either provide only one variogram to fit an isotropic model, "
            + "or directional ones for all main axes to fit anisotropy."
        )
    if is_dir_vario and model.latlon:
        raise ValueError(
            "CovModel.fit_variogram: lat-lon models don't support anisotropy."
        )
    elif model.latlon:
        # convert to yadrenko model
        x_data = 2 * np.sin(x_data / 2)
    return x_data, y_data, is_dir_vario


def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario):
    if weights is not None:
        if callable(weights):
            weights = 1.0 / weights(x_data)
        elif isinstance(weights, str) and weights == "inv":
            weights = 1.0 + x_data
        elif is_dir_vario:
            if weights.size * model.dim == x_data.size:
                weights = np.tile(weights, model.dim)
            weights = 1.0 / np.array(weights).reshape(-1)
        else:
            weights = 1.0 / np.array(weights).reshape(-1)
        curve_fit_kwargs["sigma"] = weights
        curve_fit_kwargs["absolute_sigma"] = True


def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis):
    """Create initial guess and bounds for fitting."""
    low_bounds = []
    top_bounds = []
    init_guess_list = []
    for par in DEFAULT_PARA:
        if para[par]:
            low_bounds.append(model.arg_bounds[par][0])
            if par == "var" and constrain_sill:  # var <= sill in this case
                top_bounds.append(sill)
            else:
                top_bounds.append(model.arg_bounds[par][1])
            init_guess_list.append(
                _init_guess(
                    bounds=[low_bounds[-1], top_bounds[-1]],
                    current=getattr(model, par),
                    default=model.rescale if par == "len_scale" else 1.0,
                    typ=init_guess,
                    para_name=par,
                )
            )
    for opt in model.opt_arg:
        if para[opt]:
            low_bounds.append(model.arg_bounds[opt][0])
            top_bounds.append(model.arg_bounds[opt][1])
            init_guess_list.append(
                _init_guess(
                    bounds=[low_bounds[-1], top_bounds[-1]],
                    current=getattr(model, opt),
                    default=model.default_opt_arg()[opt],
                    typ=init_guess,
                    para_name=opt,
                )
            )
    if anis:
        low_bounds += [model.anis_bounds[0]] * (model.dim - 1)
        top_bounds += [model.anis_bounds[1]] * (model.dim - 1)
        if init_guess == "default":
            def_arg = default_arg_from_bounds(model.anis_bounds)
            init_guess_list += [def_arg] * (model.dim - 1)
        elif init_guess == "current":
            init_guess_list += list(model.anis)
        else:
            raise ValueError(
                "CovModel.fit: unknown init_guess: '{}'".format(init_guess)
            )

    return (low_bounds, top_bounds), init_guess_list


def _init_guess(bounds, current, default, typ, para_name):
    """Proper determination of initial guess."""
    if typ == "default":
        if bounds[0] < default < bounds[1]:
            return default
        return default_arg_from_bounds(bounds)
    if typ == "current":
        return current
    raise ValueError("CovModel.fit: unknown init_guess: '{}'".format(typ))


def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario):
    """Create the curve for scipys curve_fit."""
    var_save = model.var

    # we need arg1, otherwise curve_fit throws an error (bug?!)
    def curve(x, arg1, *args):
        """Adapted Variogram function."""
        args = (arg1,) + args
        para_skip = 0
        opt_skip = 0
        if para["var"]:
            var_tmp = args[para_skip]
            if constrain_sill:
                nugget_tmp = sill - var_tmp
                # punishment, if resulting nugget out of range for fixed sill
                if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0:
                    return np.full_like(x, np.inf)
                # nugget estimation deselected in this case
                model.nugget = nugget_tmp
            para_skip += 1
        if para["len_scale"]:
            model.len_scale = args[para_skip]
            para_skip += 1
        if para["nugget"]:
            model.nugget = args[para_skip]
            para_skip += 1
        for opt in model.opt_arg:
            if para[opt]:
                setattr(model, opt, args[para_skip + opt_skip])
                opt_skip += 1
        # set var at last because of var_factor (other parameter needed)
        if para["var"]:
            model.var = var_tmp
        # needs to be reset for TPL models when len_scale was changed
        else:
            model.var = var_save
        if is_dir_vario:
            if anis:
                model.anis = args[1 - model.dim :]
            xs = x[: x.size // model.dim]
            out = np.array([], dtype=np.double)
            for i in range(model.dim):
                out = np.concatenate((out, model.vario_axis(xs, axis=i)))
            return out
        return model.variogram(x)

    return curve


def _post_fitting(model, para, popt, anis, is_dir_vario):
    """Postprocess fitting results and application to model."""
    fit_para = {}
    para_skip = 0
    opt_skip = 0
    for par in DEFAULT_PARA:
        if para[par]:
            if par == "var":  # set variance last
                var_tmp = popt[para_skip]
            else:
                setattr(model, par, popt[para_skip])
            fit_para[par] = popt[para_skip]
            para_skip += 1
        else:
            fit_para[par] = getattr(model, par)
    for opt in model.opt_arg:
        if para[opt]:
            setattr(model, opt, popt[para_skip + opt_skip])
            fit_para[opt] = popt[para_skip + opt_skip]
            opt_skip += 1
        else:
            fit_para[opt] = getattr(model, opt)
    if is_dir_vario:
        if anis:
            model.anis = popt[1 - model.dim :]
        fit_para["anis"] = model.anis
    # set var at last because of var_factor (other parameter needed)
    if para["var"]:
        model.var = var_tmp
    return fit_para


def _r2_score(model, x_data, y_data, is_dir_vario):
    """Calculate the R2 score."""
    if is_dir_vario:
        xs = x_data[: x_data.size // model.dim]
        vario = np.array([], dtype=np.double)
        for i in range(model.dim):
            vario = np.concatenate((vario, model.vario_axis(xs, axis=i)))
    else:
        vario = model.variogram(x_data)
    residuals = y_data - vario
    ss_res = np.sum(residuals ** 2)
    ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)
    return 1.0 - (ss_res / ss_tot)


def logistic_weights(p=0.1, mean=0.7):  # pragma: no cover
    """
    Return a logistic weights function.

    Parameters
    ----------
    p : :class:`float`, optional
        Parameter for the growth rate.
        Within this percentage of the data range, the function will
        be in the upper resp. lower percentile p. The default is 0.1.
    mean : :class:`float`, optional
        Percentage of the data range, where this function has its
        sigmoid's midpoint. The default is 0.7.

    Returns
    -------
    callable
        Weighting function.
    """
    # define the callable weights function
    def func(x_data):
        """Callable function for the weights."""
        x_range = np.amax(x_data) - np.amin(x_data)
        # logit function for growth rate
        growth = np.log(p / (1 - p)) / (p * x_range)
        x_mean = mean * x_range + np.amin(x_data)
        return 1.0 / (1.0 + np.exp(growth * (x_mean - x_data)))

    return func

1	# -- coding: utf-8 --
2	"""	15✔
3	GStools subpackage providing tools for the covariance-model.
4
5	.. currentmodule:: gstools.covmodel.fit
6
7	The following classes and functions are provided
8
9	.. autosummary::
10	fit_variogram
11	"""
12
13	# pylint: disable=C0103
14	import numpy as np	15✔
15	from scipy.optimize import curve_fit	15✔
16	from gstools.covmodel.tools import check_arg_in_bounds, default_arg_from_bounds	15✔
17	from gstools.tools.geometric import set_anis	15✔
18
19
20	__all__ = ["fit_variogram"]	15✔
21
22
23	DEFAULT_PARA = ["var", "len_scale", "nugget"]	15✔
24
25
26	def fit_variogram(	15✔
27	model,
28	x_data,
29	y_data,
30	anis=True,
31	sill=None,
32	init_guess="default",
33	weights=None,
34	method="trf",
35	loss="soft_l1",
36	max_eval=None,
37	return_r2=False,
38	curve_fit_kwargs=None,
39	**para_select
40	):
41	"""
42	Fitting a variogram-model to an empirical variogram.
43
44	Parameters
45	----------
46	model : :any:`CovModel`
47	Covariance Model to fit.
48	x_data : :class:`numpy.ndarray`
49	The bin-centers of the empirical variogram.
50	y_data : :class:`numpy.ndarray`
51	The messured variogram
52	If multiple are given, they are interpreted as the directional
53	variograms along the main axis of the associated rotated
54	coordinate system.
55	Anisotropy ratios will be estimated in that case.
56	anis : :class:`bool`, optional
57	In case of a directional variogram, you can control anisotropy
58	by this argument. Deselect the parameter from fitting, by setting
59	it "False".
60	You could also pass a fixed value to be set in the model.
61	Then the anisotropy ratios won't be altered during fitting.
62	Default: True
63	sill : :class:`float` or :class:`bool` or :any:`None`, optional
64	Here you can provide a fixed sill for the variogram.
65	It needs to be in a fitting range for the var and nugget bounds.
66	If variance or nugget are not selected for estimation,
67	the nugget will be recalculated to fulfill:
68
69	* sill = var + nugget
70	* if the variance is bigger than the sill,
71	nugget will bet set to its lower bound
72	and the variance will be set to the fitting partial sill.
73
74	If variance is deselected, it needs to be less than the sill,
75	otherwise a ValueError comes up. Same for nugget.
76	If sill=False, it will be deslected from estimation
77	and set to the current sill of the model.
78	Then, the procedure above is applied.
79	Default: None
80	init_guess : :class:`str`, optional
81	Initial guess for the estimation. Either:
82
83	* "default": using the default values of the covariance model
84	* "current": using the current values of the covariance model
85
86	Default: "default"
87	weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`optional
88	Weights applied to each point in the estimation. Either:
89
90	* 'inv': inverse distance ``1 / (x_data + 1)``
91	* list: weights given per bin
92	* callable: function applied to x_data
93
94	If callable, it must take a 1-d ndarray. Then ``weights = f(x_data)``.
95	Default: None
96	method : {'trf', 'dogbox'}, optional
97	Algorithm to perform minimization.
98
99	* 'trf' : Trust Region Reflective algorithm, particularly suitable
100	for large sparse problems with bounds. Generally robust method.
101	* 'dogbox' : dogleg algorithm with rectangular trust regions,
102	typical use case is small problems with bounds. Not recommended
103	for problems with rank-deficient Jacobian.
104
105	Default: 'trf'
106	loss : :class:`str` or :class:`callable`, optional
107	Determines the loss function in scipys curve_fit.
108	The following keyword values are allowed:
109
110	* 'linear' (default) : ``rho(z) = z``. Gives a standard
111	least-squares problem.
112	* 'soft_l1' : ``rho(z) = 2 * ((1 + z)**0.5 - 1)``. The smooth
113	approximation of l1 (absolute value) loss. Usually a good
114	choice for robust least squares.
115	* 'huber' : ``rho(z) = z if z <= 1 else 2z*0.5 - 1``. Works
116	similarly to 'soft_l1'.
117	* 'cauchy' : ``rho(z) = ln(1 + z)``. Severely weakens outliers
118	influence, but may cause difficulties in optimization process.
119	* 'arctan' : ``rho(z) = arctan(z)``. Limits a maximum loss on
120	a single residual, has properties similar to 'cauchy'.
121
122	If callable, it must take a 1-d ndarray ``z=f**2`` and return an
123	array_like with shape (3, m) where row 0 contains function values,
124	row 1 contains first derivatives and row 2 contains second
125	derivatives. Default: 'soft_l1'
126	max_eval : :class:`int` or :any:`None`, optional
127	Maximum number of function evaluations before the termination.
128	If None (default), the value is chosen automatically: 100 * n.
129	return_r2 : :class:`bool`, optional
130	Whether to return the r2 score of the estimation.
131	Default: False
132	curve_fit_kwargs : :class:`dict`, optional
133	Other keyword arguments passed to scipys curve_fit. Default: None
134	**para_select
135	You can deselect parameters from fitting, by setting
136	them "False" using their names as keywords.
137	You could also pass fixed values for each parameter.
138	Then these values will be applied and the involved parameters wont
139	be fitted.
140	By default, all parameters are fitted.
141
142	Returns
143	-------
144	fit_para : :class:`dict`
145	Dictonary with the fitted parameter values
146	pcov : :class:`numpy.ndarray`
147	The estimated covariance of `popt` from
148	:any:`scipy.optimize.curve_fit`.
149	To compute one standard deviation errors
150	on the parameters use ``perr = np.sqrt(np.diag(pcov))``.
151	r2_score : :class:`float`, optional
152	r2 score of the curve fitting results. Only if return_r2 is True.
153
154	Notes
155	-----
156	You can set the bounds for each parameter by accessing
157	:any:`CovModel.set_arg_bounds`.
158
159	The fitted parameters will be instantly set in the model.
160	"""
161	# preprocess selected parameters
162	para, sill, constrain_sill, anis = _pre_para(	15✔
163	model, para_select, sill, anis
164	)
165	# check curve_fit kwargs
166	curve_fit_kwargs = {} if curve_fit_kwargs is None else curve_fit_kwargs	15✔
167	# check method
168	if method not in ["trf", "dogbox"]:	15✔
169	raise ValueError("fit: method needs to be either 'trf' or 'dogbox'")	15✔
170	# prepare variogram data
171	# => concatenate directional variograms to have a 1D array for x and y
172	x_data, y_data, is_dir_vario = _check_vario(model, x_data, y_data)	15✔
173	# only fit anisotropy if a directional variogram was given
174	anis &= is_dir_vario	15✔
175	# set weights
176	_set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario)	15✔
177	# set the lower/upper boundaries for the variogram-parameters
178	bounds, init_guess_list = _init_curve_fit_para(	15✔
179	model, para, init_guess, constrain_sill, sill, anis
180	)
181	# create the fitting curve
182	curve_fit_kwargs["f"] = _get_curve(	15✔
183	model, para, constrain_sill, sill, anis, is_dir_vario
184	)
185	# set the remaining kwargs for curve_fit
186	curve_fit_kwargs["bounds"] = bounds	15✔
187	curve_fit_kwargs["p0"] = init_guess_list	15✔
188	curve_fit_kwargs["xdata"] = x_data	15✔
189	curve_fit_kwargs["ydata"] = y_data	15✔
190	curve_fit_kwargs["loss"] = loss	15✔
191	curve_fit_kwargs["max_nfev"] = max_eval	15✔
192	curve_fit_kwargs["method"] = method	15✔
193	# fit the variogram
194	popt, pcov = curve_fit(**curve_fit_kwargs)	15✔
195	# convert the results
196	fit_para = _post_fitting(model, para, popt, anis, is_dir_vario)	15✔
197	# calculate the r2 score if wanted
198	if return_r2:	15✔
199	return fit_para, pcov, _r2_score(model, x_data, y_data, is_dir_vario)	15✔
200	return fit_para, pcov	15✔
201
202
203	def _pre_para(model, para_select, sill, anis):	15✔
204	"""Preprocess selected parameters."""
205	var_last = False	15✔
206	for par in para_select:	15✔
207	if par not in model.arg_bounds:	15✔
208	raise ValueError(	15✔
209	"fit: unknow parameter in selection: {}".format(par)
210	)
211	if not isinstance(para_select[par], bool):	15✔
212	if par == "var":	15✔
213	var_last = True	15✔
214	var_tmp = float(para_select[par])	15✔
215	else:
216	setattr(model, par, float(para_select[par]))	15✔
217	para_select[par] = False	15✔
218	# set variance last due to possible recalculations
219	if var_last:	15✔
220	model.var = var_tmp	15✔
221	# remove those that were set to True
222	para_select = {k: v for k, v in para_select.items() if not v}	15✔
223	# handling the sill
224	sill = None if (isinstance(sill, bool) and sill) else sill	15✔
225	if sill is not None:	15✔
226	sill = model.sill if isinstance(sill, bool) else float(sill)	15✔
227	constrain_sill = True	15✔
228	sill_low = model.arg_bounds["var"][0] + model.arg_bounds["nugget"][0]	15✔
229	sill_up = model.arg_bounds["var"][1] + model.arg_bounds["nugget"][1]	15✔
230	if not sill_low <= sill <= sill_up:	15✔
231	raise ValueError("fit: sill out of bounds.")	15✔
232	if "var" in para_select and "nugget" in para_select:	15✔
233	if model.var > sill:	15✔
234	model.nugget = model.arg_bounds["nugget"][0]	15✔
235	model.var = sill - model.nugget	15✔
236	else:
237	model.nugget = sill - model.var	15✔
238	elif "var" in para_select:	15✔
239	if model.var > sill:	15✔
240	raise ValueError(	15✔
241	"fit: if sill is fixed and variance deselected, "
242	+ "the set variance should be less than the given sill."
243	)
244	para_select["nugget"] = False	15✔
245	model.nugget = sill - model.var	15✔
246	elif "nugget" in para_select:	15✔
247	if model.nugget > sill:	15✔
248	raise ValueError(	15✔
249	"fit: if sill is fixed and nugget deselected, "
250	+ "the set nugget should be less than the given sill."
251	)
252	para_select["var"] = False	15✔
253	model.var = sill - model.nugget	15✔
254	else:
255	# deselect the nugget, to recalculate it accordingly
256	# nugget = sill - var
257	para_select["nugget"] = False	15✔
258	else:
259	constrain_sill = False	15✔
260	# select all parameters to be fitted
261	para = {par: True for par in DEFAULT_PARA}	15✔
262	para.update({opt: True for opt in model.opt_arg})	15✔
263	# now deselect unwanted parameters
264	para.update(para_select)	15✔
265	# check if anisotropy should be fitted or set
266	if not isinstance(anis, bool):	15✔
267	model.anis = anis	15✔
268	anis = False	15✔
269	return para, sill, constrain_sill, anis	15✔
270
271
272	def _check_vario(model, x_data, y_data):	15✔
273	# prepare variogram data
274	x_data = np.array(x_data).reshape(-1)	15✔
275	y_data = np.array(y_data).reshape(-1)	15✔
276	# if multiple variograms are given, they will be interpreted
277	# as directional variograms along the main rotated axes of the model
278	is_dir_vario = False	15✔
279	if model.dim > 1 and x_data.size * model.dim == y_data.size:	15✔
280	is_dir_vario = True	15✔
281	# concatenate multiple variograms
282	x_data = np.tile(x_data, model.dim)	15✔
283	elif x_data.size != y_data.size:	15✔
284	raise ValueError(	15✔
285	"CovModel.fit_variogram: Wrong number of empirical variograms! "
286	+ "Either provide only one variogram to fit an isotropic model, "
287	+ "or directional ones for all main axes to fit anisotropy."
288	)
289	if is_dir_vario and model.latlon:	15✔
290	raise ValueError(	×
291	"CovModel.fit_variogram: lat-lon models don't support anisotropy."
292	)
293	elif model.latlon:	15✔
294	# convert to yadrenko model
295	x_data = 2 * np.sin(x_data / 2)	15✔
296	return x_data, y_data, is_dir_vario	15✔
297
298
299	def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario):	15✔
300	if weights is not None:	15✔
301	if callable(weights):	15✔
302	weights = 1.0 / weights(x_data)	15✔
303	elif isinstance(weights, str) and weights == "inv":	15✔
304	weights = 1.0 + x_data	15✔
305	elif is_dir_vario:	15✔
306	if weights.size * model.dim == x_data.size:	15✔
307	weights = np.tile(weights, model.dim)	15✔
308	weights = 1.0 / np.array(weights).reshape(-1)	15✔
309	else:
310	weights = 1.0 / np.array(weights).reshape(-1)	×
311	curve_fit_kwargs["sigma"] = weights	15✔
312	curve_fit_kwargs["absolute_sigma"] = True	15✔
313
314
315	def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis):	15✔
316	"""Create initial guess and bounds for fitting."""
317	low_bounds = []	15✔
318	top_bounds = []	15✔
319	init_guess_list = []	15✔
320	for par in DEFAULT_PARA:	15✔
321	if para[par]:	15✔
322	low_bounds.append(model.arg_bounds[par][0])	15✔
323	if par == "var" and constrain_sill: # var <= sill in this case	15✔
324	top_bounds.append(sill)	15✔
325	else:
326	top_bounds.append(model.arg_bounds[par][1])	15✔
327	init_guess_list.append(	15✔
328	_init_guess(
329	bounds=[low_bounds[-1], top_bounds[-1]],
330	current=getattr(model, par),
331	default=model.rescale if par == "len_scale" else 1.0,
332	typ=init_guess,
333	para_name=par,
334	)
335	)
336	for opt in model.opt_arg:	15✔
337	if para[opt]:	15✔
338	low_bounds.append(model.arg_bounds[opt][0])	15✔
339	top_bounds.append(model.arg_bounds[opt][1])	15✔
340	init_guess_list.append(	15✔
341	_init_guess(
342	bounds=[low_bounds[-1], top_bounds[-1]],
343	current=getattr(model, opt),
344	default=model.default_opt_arg()[opt],
345	typ=init_guess,
346	para_name=opt,
347	)
348	)
349	if anis:	15✔
350	low_bounds += [model.anis_bounds[0]] * (model.dim - 1)	15✔
351	top_bounds += [model.anis_bounds[1]] * (model.dim - 1)	15✔
352	if init_guess == "default":	15✔
353	def_arg = default_arg_from_bounds(model.anis_bounds)	15✔
354	init_guess_list += [def_arg] * (model.dim - 1)	15✔
355	elif init_guess == "current":	×
356	init_guess_list += list(model.anis)	×
357	else:
358	raise ValueError(	×
359	"CovModel.fit: unknown init_guess: '{}'".format(init_guess)
360	)
361
362	return (low_bounds, top_bounds), init_guess_list	15✔
363
364
365	def _init_guess(bounds, current, default, typ, para_name):	15✔
366	"""Proper determination of initial guess."""
367	if typ == "default":	15✔
368	if bounds[0] < default < bounds[1]:	15✔
369	return default	15✔
370	return default_arg_from_bounds(bounds)	15✔
371	if typ == "current":	15✔
372	return current	15✔
373	raise ValueError("CovModel.fit: unknown init_guess: '{}'".format(typ))	15✔
374
375
376	def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario):	15✔
377	"""Create the curve for scipys curve_fit."""
378	var_save = model.var	15✔
379
380	# we need arg1, otherwise curve_fit throws an error (bug?!)
381	def curve(x, arg1, *args):	15✔
382	"""Adapted Variogram function."""
383	args = (arg1,) + args	15✔
384	para_skip = 0	15✔
385	opt_skip = 0	15✔
386	if para["var"]:	15✔
387	var_tmp = args[para_skip]	15✔
388	if constrain_sill:	15✔
389	nugget_tmp = sill - var_tmp	15✔
390	# punishment, if resulting nugget out of range for fixed sill
391	if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0:	15✔
392	return np.full_like(x, np.inf)	15✔
393	# nugget estimation deselected in this case
394	model.nugget = nugget_tmp	15✔
395	para_skip += 1	15✔
396	if para["len_scale"]:	15✔
397	model.len_scale = args[para_skip]	15✔
398	para_skip += 1	15✔
399	if para["nugget"]:	15✔
400	model.nugget = args[para_skip]	15✔
401	para_skip += 1	15✔
402	for opt in model.opt_arg:	15✔
403	if para[opt]:	15✔
404	setattr(model, opt, args[para_skip + opt_skip])	15✔
405	opt_skip += 1	15✔
406	# set var at last because of var_factor (other parameter needed)
407	if para["var"]:	15✔
408	model.var = var_tmp	15✔
409	# needs to be reset for TPL models when len_scale was changed
410	else:
411	model.var = var_save	15✔
412	if is_dir_vario:	15✔
413	if anis:	15✔
414	model.anis = args[1 - model.dim :]	15✔
415	xs = x[: x.size // model.dim]	15✔
416	out = np.array([], dtype=np.double)	15✔
417	for i in range(model.dim):	15✔
418	out = np.concatenate((out, model.vario_axis(xs, axis=i)))	15✔
419	return out	15✔
420	return model.variogram(x)	15✔
421
422	return curve	15✔
423
424
425	def _post_fitting(model, para, popt, anis, is_dir_vario):	15✔
426	"""Postprocess fitting results and application to model."""
427	fit_para = {}	15✔
428	para_skip = 0	15✔
429	opt_skip = 0	15✔
430	for par in DEFAULT_PARA:	15✔
431	if para[par]:	15✔
432	if par == "var": # set variance last	15✔
433	var_tmp = popt[para_skip]	15✔
434	else:
435	setattr(model, par, popt[para_skip])	15✔
436	fit_para[par] = popt[para_skip]	15✔
437	para_skip += 1	15✔
438	else:
439	fit_para[par] = getattr(model, par)	15✔
440	for opt in model.opt_arg:	15✔
441	if para[opt]:	15✔
442	setattr(model, opt, popt[para_skip + opt_skip])	15✔
443	fit_para[opt] = popt[para_skip + opt_skip]	15✔
444	opt_skip += 1	15✔
445	else:
446	fit_para[opt] = getattr(model, opt)	15✔
447	if is_dir_vario:	15✔
448	if anis:	15✔
449	model.anis = popt[1 - model.dim :]	15✔
450	fit_para["anis"] = model.anis	15✔
451	# set var at last because of var_factor (other parameter needed)
452	if para["var"]:	15✔
453	model.var = var_tmp	15✔
454	return fit_para	15✔
455
456
457	def _r2_score(model, x_data, y_data, is_dir_vario):	15✔
458	"""Calculate the R2 score."""
459	if is_dir_vario:	15✔
460	xs = x_data[: x_data.size // model.dim]	15✔
461	vario = np.array([], dtype=np.double)	15✔
462	for i in range(model.dim):	15✔
463	vario = np.concatenate((vario, model.vario_axis(xs, axis=i)))	15✔
464	else:
465	vario = model.variogram(x_data)	15✔
466	residuals = y_data - vario	15✔
467	ss_res = np.sum(residuals ** 2)	15✔
468	ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)	15✔
469	return 1.0 - (ss_res / ss_tot)	15✔
470
471
472	def logistic_weights(p=0.1, mean=0.7): # pragma: no cover
473	"""
474	Return a logistic weights function.
475
476	Parameters
477	----------
478	p : :class:`float`, optional
479	Parameter for the growth rate.
480	Within this percentage of the data range, the function will
481	be in the upper resp. lower percentile p. The default is 0.1.
482	mean : :class:`float`, optional
483	Percentage of the data range, where this function has its
484	sigmoid's midpoint. The default is 0.7.
485
486	Returns
487	-------
488	callable
489	Weighting function.
490	"""
491	# define the callable weights function
492	def func(x_data):
493	"""Callable function for the weights."""
494	x_range = np.amax(x_data) - np.amin(x_data)
495	# logit function for growth rate
496	growth = np.log(p / (1 - p)) / (p * x_range)
497	x_mean = mean * x_range + np.amin(x_data)
498	return 1.0 / (1.0 + np.exp(growth * (x_mean - x_data)))
499
500	return func

GeoStat-Framework / GSTools / 476491407

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous