Variogram.py

"""
Variogram class
"""
import copy
import os
import warnings

import numpy as np
from pandas import DataFrame
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.spatial.distance import pdist, squareform
from sklearn.isotonic import IsotonicRegression

from skgstat import estimators, models, binning


class Variogram(object):
    """Variogram Class

    Calculates a variogram of the separating distances in the given
    coordinates and relates them to one of the semi-variance measures of
    the given dependent values.

    """
    def __init__(self,
                 coordinates=None,
                 values=None,
                 estimator='matheron',
                 model='spherical',
                 dist_func='euclidean',
                 bin_func='even',
                 normalize=False,
                 fit_method='trf',
                 fit_sigma=None,
                 use_nugget=False,
                 maxlag=None,
                 n_lags=10,
                 verbose=False
                 ):
        r"""Variogram Class

        Note: The directional variogram estimation is not re-implemented yet.
        Therefore the parameters is-directional, azimuth and tolerance will
        be ignored at the moment and can be subject to changes.

        Parameters
        ----------
        coordinates : numpy.ndarray
            Array of shape (m, n). Will be used as m observation points of
            n-dimensions. This variogram can be calculated on 1 - n
            dimensional coordinates. In case a 1-dimensional array is passed,
            a second array of same length containing only zeros will be
            stacked to the passed one.
        values : numpy.ndarray
            Array of values observed at the given coordinates. The length of
            the values array has to match the m dimension of the coordinates
            array. Will be used to calculate the dependent variable of the
            variogram.
        estimator : str, callable
            String identifying the semi-variance estimator to be used.
            Defaults to the Matheron estimator. Possible values are:

              * matheron        [Matheron, default]
              * cressie         [Cressie-Hawkins]
              * dowd            [Dowd-Estimator]
              * genton          [Genton]
              * minmax          [MinMax Scaler]
              * entropy         [Shannon Entropy]

            If a callable is passed, it has to accept an array of absoulte
            differences, aligned to the 1D distance matrix (flattened upper
            triangle) and return a scalar, that converges towards small
            values for similarity (high covariance).
        model : str
            String identifying the theoretical variogram function to be used
            to describe the experimental variogram. Can be one of:

              * spherical       [Spherical, default]
              * exponential     [Exponential]
              * gaussian        [Gaussian]
              * cubic           [Cubic]
              * stable          [Stable model]
              * matern          [Matérn model]
              * nugget          [nugget effect variogram]

        dist_func : str
            String identifying the distance function. Defaults to
            'euclidean'. Can be any metric accepted by
            scipy.spatial.distance.pdist. Additional parameters are not (yet)
            passed through to pdist. These are accepted by pdist for some of
            the metrics. In these cases the default values are used.
        bin_func : str
            String identifying the binning function used to find lag class
            edges. At the moment there are two possible values: 'even'
            (default) or 'uniform'. Even will find n_lags bins of same width
            in the interval [0,maxlag[. 'uniform' will identfy n_lags bins on
            the same interval, but with varying edges so that all bins count
            the same amount of observations.
        normalize : bool
            Defaults to False. If True, the independent and dependent
            variable will be normalized to the range [0,1].
        fit_method : str
            String identifying the method to be used for fitting the
            theoretical variogram function to the experimental. More info is
            given in the Variogram.fit docs. Can be one of:

                * 'lm': Levenberg-Marquardt algorithm for unconstrained
                  problems. This is the faster algorithm, yet is the fitting of
                  a variogram not unconstrianed.
                * 'trf': Trust Region Reflective function for non-linear
                  constrained problems. The class will set the boundaries
                  itself. This is the default function.

        fit_sigma : numpy.ndarray, str
            Defaults to None. The sigma is used as measure of uncertainty
            during variogram fit. If fit_sigma is an array, it has to hold
            n_lags elements, giving the uncertainty for all lags classes. If
            fit_sigma is None (default), it will give no weight to any lag.
            Higher values indicate higher uncertainty and will lower the
            influcence of the corresponding lag class for the fit.
            If fit_sigma is a string, a pre-defined function of separating
            distance will be used to fill the array. Can be one of:

                * 'linear': Linear loss with distance. Small bins will have
                  higher impact.
                * 'exp': The weights decrease by a e-function of distance
                * 'sqrt': The weights decrease by the squareroot of distance
                * 'sq': The weights decrease by the squared distance.

            More info is given in the Variogram.fit_sigma documentation.
        use_nugget : bool
            Defaults to False. If True, a nugget effet will be added to all
            Variogram.models as a third (or fourth) fitting parameter. A
            nugget is essentially the y-axis interception of the theoretical
            variogram function.
        maxlag : float, str
            Can specify the maximum lag distance directly by giving a value
            larger than 1. The binning function will not find any lag class
            with an edge larger than maxlag. If 0 < maxlag < 1, then maxlag
            is relative and maxlag * max(Variogram.distance) will be used.
            In case maxlag is a string it has to be one of 'median', 'mean'.
            Then the median or mean of all Variogram.distance will be used.
            Note maxlag=0.5 will use half the maximum separating distance,
            this is not the same as 'median', which is the median of all
            separating distances
        n_lags : int
            Specify the number of lag classes to be defined by the binning
            function.
        verbose : bool
            Set the Verbosity of the class. Not Implemented yet.

        """
        # Set coordinates
        self._X = np.asarray(coordinates)

        # pairwise differences
        self._diff = None

        # set verbosity
        self.verbose = verbose

        # set values
        self._values = None
        # calc_diff = False here, because it will be calculated by fit() later
        self.set_values(values=values, calc_diff=False)

        # distance matrix
        self._dist = None

        # set distance calculation function
        self._dist_func = None
        self.set_dist_function(func=dist_func)

        # lags and max lag
        self._n_lags = None
        self.n_lags = n_lags
        self._maxlag = None
        self.maxlag = maxlag

        # harmonize model placeholder
        self._harmonize = False

        # estimator can be a function or a string
        self._estimator = None
        self.set_estimator(estimator_name=estimator)

        # model can be a function or a string
        self._model = None
        self.set_model(model_name=model)

        # the binning settings
        self._bin_func = None
        self._groups = None
        self._bins = None
        self.set_bin_func(bin_func=bin_func)

        # specify if the lag should be given absolute or relative to the maxlag
        self._normalized = normalize

        # set if nugget effect shall be used
        self._use_nugget = None
        self.use_nugget = use_nugget

        # set the fitting method and sigma array
        self.fit_method = fit_method
        self._fit_sigma = None
        self.fit_sigma = fit_sigma

        # set attributes to be filled during calculation
        self.cov = None
        self.cof = None

        # settings, not reachable by init (not yet)
        self._cache_experimental = False

        # do the preprocessing and fitting upon initialization
        # Note that fit() calls preprocessing
        self.fit(force=True)

    @property
    def coordinates(self):
        """Coordinates property

        Array of observation locations the variogram is build for. This
        property has no setter. If you want to change the coordinates,
        use a new Variogram instance.

        Returns
        -------
        coordinates : numpy.array

        """
        return self._X

    @property
    def values(self):
        """Values property

        Array of observations, the variogram is build for. The setter of this
        property utilizes the Variogram.set_values function for setting new
        arrays.

        Returns
        -------
        values : numpy.ndarray

        See Also
        --------
        Variogram.set_values

        """
        return self._values

    @values.setter
    def values(self, values):
        self.set_values(values=values)

    @property
    def value_matrix(self):
        """Value matrix

        Returns a matrix of pairwise differences in absolute values. The
        matrix will have the shape (m, m) with m = len(Variogram.values).
        Note that Variogram.values holds the values themselves, while the
        value_matrix consists of their pairwise differences.

        Returns
        -------
        values : numpy.matrix
            Matrix of pairwise absolute differences of the values.

        See Also
        --------
        Variogram._diff

        """
        return squareform(self._diff)

    def set_values(self, values, calc_diff=True):
        """Set new values

        Will set the passed array as new value array. This array has to be of
        same length as the first axis of the coordinates array. The Variogram
        class does only accept one dimensional arrays.
        On success all fitting parameters are deleted and the pairwise
        differences are recalculated.
        Raises :py:class:`ValueError`s on shape mismatches and a Warning

        Parameters
        ----------
        values : numpy.ndarray

        Returns
        -------
        void

        Raises
        ------
        ValueError : raised if the values array shape does not match the
            coordinates array, or more than one dimension given
        Warning : raised if all input values are the same

        See Also
        --------
        Variogram.values

        """
        # check dimensions
        if not len(values) == len(self._X):
            raise ValueError('The length of the values array has to match' +
                             'the length of coordinates')

        # use an array
        _y = np.asarray(values)
        if not _y.ndim == 1:
            raise ValueError('The values shall be a 1-D array.' +
                             'Multi-dimensional values not supported yet.')

        # check if all input values are the same
        if len(set(_y)) < 2:
            raise Warning('All input values are the same.')

        # reset fitting parameter
        self.cof, self.cov = None, None
        self._diff = None

        # set new values
        self._values = np.asarray(values)

        # recalculate the pairwise differences
        if calc_diff:
            self._calc_diff(force=True)

    @property
    def bin_func(self):
        """Binning function

        Returns an instance of the function used for binning the separating
        distances into the given amount of bins. Both functions use the same
        signature of func(distances, n, maxlag).

        The setter of this property utilizes the Variogram.set_bin_func to
        set a new function.

        Returns
        -------
        binning_function : function

        See Also
        --------
        Variogram.set_bin_func

        """
        return self._bin_func

    @bin_func.setter
    def bin_func(self, bin_func):
        self.set_bin_func(bin_func=bin_func)

    def set_bin_func(self, bin_func):
        """Set binning function

        Sets a new binning function to be used. The new binning method is set
        by a string identifying the new function to be used. Can be one of:
        ['even', 'uniform'].

        Parameters
        ----------
        bin_func : str
            Can be one of:

            * **'even'**: Use skgstat.binning.even_width_lags for using
              n_lags lags of equal width up to maxlag.
            * **'uniform'**: Use skgstat.binning.uniform_count_lags for using
              n_lags lags up to maxlag in which the pairwise differences
              follow a uniform distribution.

        Returns
        -------
        void

        See Also
        --------
        Variogram.bin_func
        skgstat.binning.uniform_count_lags
        skgstat.binning.even_width_lags

        """
        # switch the input
        if bin_func.lower() == 'even':
            self._bin_func = binning.even_width_lags
        elif bin_func.lower() == 'uniform':
            self._bin_func = binning.uniform_count_lags
        else:
            raise ValueError('%s binning method is not known' % bin_func)

        # reset groups and bins
        self._groups = None
        self._bins = None
        self.cof, self.cov = None, None

    @property
    def normalized(self):
        return self._normalized

    @normalized.setter
    def normalized(self, status):
        # set the new value
        self._normalized = status

    @property
    def bins(self):
        if self._bins is None:
            self._bins = self.bin_func(self.distance, self.n_lags, self.maxlag)

        return self._bins.copy()

    @bins.setter
    def bins(self, bins):
        # set the new bins
        self._bins = bins

        # clean the groups as they are not valid anymore
        self._groups = None
        self.cov = None
        self.cof = None

    @property
    def n_lags(self):
        """Number of lag bins

        Pass the number of lag bins to be used on
        this Variogram instance. This will reset
        the grouping index and fitting parameters

        """
        return self._n_lags

    @n_lags.setter
    def n_lags(self, n):
        # TODO: here accept strings and implement some optimum methods
        # string are not implemented yet
        if isinstance(n, str):
            raise NotImplementedError('n_lags string values not implemented')

        # n_lags is int
        elif isinstance(n, int):
            if n < 1:
                raise ValueError('n_lags has to be a positive integer')

            # set parameter
            self._n_lags = n

            # reset the bins
            self._bins = None

        # else
        else:
            raise ValueError('n_lags has to be a positive integer')

        # reset the fitting
        self.cof = None
        self.cov = None

    @property
    def estimator(self):
        return self._estimator

    @estimator.setter
    def estimator(self, value):
        self.set_estimator(estimator_name=value)

    def set_estimator(self, estimator_name):
        # reset the fitting
        self.cof, self.cov = None, None

        if isinstance(estimator_name, str):
            if estimator_name.lower() == 'matheron':
                self._estimator = estimators.matheron
            elif estimator_name.lower() == 'cressie':
                self._estimator = estimators.cressie
            elif estimator_name.lower() == 'dowd':
                self._estimator = estimators.dowd
            elif estimator_name.lower() == 'genton':
                self._estimator = estimators.genton
            elif estimator_name.lower() == 'minmax':
                self._estimator = estimators.minmax
            elif estimator_name.lower() == 'percentile':
                self._estimator = estimators.percentile
            elif estimator_name.lower() == 'entropy':
                self._estimator = estimators.entropy
            else:
                raise ValueError(
                    ('Variogram estimator %s is not understood, please ' +
                    'provide the function.') % estimator_name
                )
        elif callable(estimator_name):
            self._estimator = estimator_name
        else:
            raise ValueError('The estimator has to be a string or callable.')

    @property
    def model(self):
        return self._model

    @model.setter
    def model(self, value):
        self.set_model(model_name=value)

    def set_model(self, model_name):
        """
        Set model as the new theoretical variogram function.

        """
        # reset the fitting
        self.cof, self.cov = None, None

        if isinstance(model_name, str):
            # at first reset harmonize
            self._harmonize = False
            if model_name.lower() == 'spherical':
                self._model = models.spherical
            elif model_name.lower() == 'exponential':
                self._model = models.exponential
            elif model_name.lower() == 'gaussian':
                self._model = models.gaussian
            elif model_name.lower() == 'cubic':
                self._model = models.cubic
            elif model_name.lower() == 'stable':
                self._model = models.stable
            elif model_name.lower() == 'matern':
                self._model = models.matern
            elif model_name.lower() == 'harmonize':
                self._harmonize = True
                self._model = self._build_harmonized_model()
            else:
                raise ValueError(
                    'The theoretical Variogram function %s is not understood, please provide the function' % model_name)
        else:
            self._model = model_name

    def _build_harmonized_model(self):
        x = self.bins
        y = self.experimental

        _x = x[~np.isnan(y)]
        _y = y[~np.isnan(y)]
        regr = IsotonicRegression(increasing=True).fit(_x, _y)

        # create the model function
        def harmonize(x):
            """Monotonized Variogram

            Return the isotonic harmonized experimental variogram.
            This means, the experimental variogram is monotonic after harmonization.

            The harmonization is done using following Hinterding (2003) using 
            the PAVA algorithm (Barlow and Bartholomew, 1972).

            Returns
            -------
            gamma : numpy.ndarray
                monotonized experimental variogram
            
            References
            ----------
            Barlow, R., D. Bartholomew, et al. (1972): Statistical Interference Under Order Restrictions.
                John Wiley and Sons, New York.
            Hiterding, A. (2003): Entwicklung hybrider Interpolationsverfahren für den automatisierten Betrieb am
                Beispiel meteorologischer Größen. Dissertation, Institut für Geoinformatik, Westphälische
                Wilhelms-Universität Münster, IfGIprints, Münster. ISBN: 3-936616-12-4

            """

            if isinstance(x, (list, tuple, np.ndarray)):
                return regr.transform(x)
            else:
                return regr.transform([x])

        return harmonize

    @property
    def use_nugget(self):
        return self._use_nugget

    @use_nugget.setter
    def use_nugget(self, nugget):
        if not isinstance(nugget, bool):
            raise ValueError('use_nugget has to be of type bool.')

        # set new value
        self._use_nugget = nugget

    @property
    def dist_function(self):
        return self._dist_func

    def _dist_func_wrapper(self, x):
        if callable(self._dist_func):
            return self._dist_func(x)
        else:
            return pdist(X=x, metric=self._dist_func)
    
    @dist_function.setter
    def dist_function(self, func):
        self.set_dist_function(func=func)

    def set_dist_function(self, func):
        """Set distance function

        Set the function used for distance calculation. func can either be a
        callable or a string. The ranked distance function is not implemented
        yet. strings will be forwarded to the scipy.spatial.distance.pdist
        function as the metric argument.
        If func is a callable, it has to return the upper triangle of the
        distance matrix as a flat array (Like the pdist function).

        Parameters
        ----------
        func : string, callable

        Returns
        -------
        numpy.array

        """
        # reset the distances and fitting
        self._dist = None
        self.cof, self.cov = None, None

        if isinstance(func, str):
            if func.lower() == 'rank':
                raise NotImplementedError
            else:
                # if not ranks, it has to be a scipy metric
                self._dist_func = func

        elif callable(func):
            self._dist_func = func
        else:
            raise ValueError('Input not supported. Pass a string or callable.')

        # re-calculate distances
        self._calc_distances()

    @property
    def distance(self):
        if self._dist is None:
            self._calc_distances()
        return self._dist

    @distance.setter
    def distance(self, dist_array):
        self._dist = dist_array

    @property
    def distance_matrix(self):
        return squareform(self.distance)

    @property
    def maxlag(self):
        return self._maxlag

    @maxlag.setter
    def maxlag(self, value):
        # reset fitting
        self.cof, self.cov = None, None

        # remove bins
        self._bins = None
        self._groups = None

        # set new maxlag
        if value is None:
            self._maxlag = None
        elif isinstance(value, str):
            if value == 'median':
                self._maxlag = np.median(self.distance)
            elif value == 'mean':
                self._maxlag = np.mean(self.distance)
        elif value < 1:
            self._maxlag = value * np.max(self.distance)
        else:
            self._maxlag = value

    @property
    def fit_sigma(self):
        r"""Fitting Uncertainty

        Set or calculate an array of observation uncertainties aligned to the
        Variogram.bins. These will be used to weight the observations in the
        cost function, which divides the residuals by their uncertainty.

        When setting fit_sigma, the array of uncertainties itself can be
        given, or one of the strings: ['linear', 'exp', 'sqrt', 'sq']. The
        parameters described below refer to the setter of this property.

        Parameters
        ----------
        sigma : string, array
            Sigma can either be an array of discrete uncertainty values,
            which have to align to the Variogram.bins, or of type string.
            Then, the weights for fitting are calculated as a function of
            (lag) distance.

              * **sigma='linear'**: The residuals get weighted by the lag
                distance normalized to the maximum lag distance, denoted as
                :math:`w_n`
              * **sigma='exp'**: The residuals get weighted by the function:
                :math:`w = e^{1 / w_n}`
              * **sigma='sqrt'**: The residuals get weighted by the function:
                :math:`w = \sqrt(w_n)`
              * **sigma='sq'**: The residuals get weighted by the function:
                :math:`w = w_n^2`

        Returns
        -------
        void

        Notes
        -----
        The cost function is defined as:

        .. math::
            chisq = \sum {\frac{r}{\sigma}}^2

        where r are the residuals between the experimental variogram and the
        modeled values for the same lag. Following this function,
        small values will increase the influence of that residual, while a very
        large sigma will cause the observation to be ignored.

        See Also
        --------
        scipy.optimize.curve_fit

        """
        # unweighted
        if self._fit_sigma is None:
            return None

        # discrete uncertainties
        elif isinstance(self._fit_sigma, (list, tuple, np.ndarray)):
            if len(self._fit_sigma) == len(self._bins):
                return self._fit_sigma
            else:
                raise AttributeError('fit_sigma and bins need the same length.')

        # linear function of distance
        elif self._fit_sigma == 'linear':
            return self.bins / np.max(self.bins)

        # e function of distance
        elif self._fit_sigma == 'exp':
            return 1. / np.exp(1. / (self.bins / np.max(self.bins)))

        # sqrt function of distance
        elif self._fit_sigma == 'sqrt':
            return np.sqrt(self.bins / np.max(self.bins))

        # squared function of distance
        elif self._fit_sigma == 'sq':
            return (self.bins / np.max(self.bins)) ** 2
        else:
            raise ValueError("fit_sigma is not understood. It has to be an " +
                             "array or one of ['linear', 'exp', 'sqrt', 'sq'].")

    @fit_sigma.setter
    def fit_sigma(self, sigma):
        self._fit_sigma = sigma

        # remove fitting parameters
        self.cof = None
        self.cov = None

    def lag_groups(self):
        """Lag class groups

        Retuns a mask array with as many elements as self._diff has,
        identifying the lag class group for each pairwise difference. Can be
        used to extract all pairwise values within the same lag bin.

        Returns
        -------
        numpy.ndarray

        See Also
        --------
        Variogram.lag_classes

        """
        if self._groups is None:
            self._calc_groups()

        return self._groups

    def lag_classes(self):
        """Iterate over the lag classes

        Generates an iterator over all lag classes. Can be zipped with
        Variogram.bins to identify the lag.

        Returns
        -------
        iterable

        """
        # yield all groups
        for i in np.unique(self.lag_groups()):
            if i < 0:
                continue
            else:
                yield self._diff[np.where(self.lag_groups() == i)]

    def preprocessing(self, force=False):
        """Preprocessing function

        Prepares all input data for the fit and transform functions. Namely,
        the distances are calculated and the value differences. Then the
        binning is set up and bin edges are calculated. If any of the listed
        subsets are already prepared, their processing is skipped. This
        behaviour can be changed by the force parameter. This will cause a
        clean preprocessing.

        Parameters
        ----------
        force : bool
            If set to True, all preprocessing data sets will be deleted. Use
            it in case you need a clean preprocessing.

        Returns
        -------
        void

        """
        # call the _calc functions
        self._calc_distances(force=force)
        self._calc_diff(force=force)
        self._calc_groups(force=force)

    def fit(self, force=False, method=None, sigma=None, **kwargs):
        """Fit the variogram

        The fit function will fit the theoretical variogram function to the
        experimental. The preprocessed distance matrix, pairwise differences
        and binning will not be recalculated, if already done. This could be
        forced by setting the force parameter to true.

        In case you call fit function directly, with method or sigma,
        the parameters set on Variogram object instantiation will get
        overwritten. All other keyword arguments will be passed to
        scipy.optimize.curve_fit function.

        Parameters
        ----------
        force : bool
            If set to True, a clean preprocessing of the distance matrix,
            pairwise differences and the binning will be forced. Default is
            False.
        method : string
            A string identifying one of the implemented fitting procedures.
            Can be one of ['lm', 'trf']:

              * lm: Levenberg-Marquardt algorithms implemented in
                scipy.optimize.leastsq function.
              * trf: Trust Region Reflective algorithm implemented in
                scipy.optimize.least_squares(method='trf')

        sigma : string, array
            Uncertainty array for the bins. Has to have the same dimension as
            self.bins. Refer to Variogram.fit_sigma for more information.

        Returns
        -------
        void

        See Also
        --------
        scipy.optimize
        scipy.optimize.curve_fit
        scipy.optimize.leastsq
        scipy.optimize.least_squares

        """
        # TODO: the kwargs need to be preserved somehow
        # delete the last cov and cof
        self.cof = None
        self.cov = None

        # if force, force a clean preprocessing
        self.preprocessing(force=force)

        # load the data
        x = self.bins
        y = self.experimental

        # overwrite fit setting if new params are given
        if method is not None:
            self.fit_method = method
        if sigma is not None:
            self.fit_sigma = sigma

        # remove nans
        _x = x[~np.isnan(y)]
        _y = y[~np.isnan(y)]

        # handle harmonized models
        if self._harmonize:
            _x = np.linspace(0, np.max(_x), 100)
            _y = self._model(_x)

            # get the params
            s = 0.95 * np.nanmax(_y)
            with warnings.catch_warnings():
                warnings.simplefilter('ignore')
                r = np.argwhere(_y >= s)[0][0]
            n = _y[0] if self.use_nugget else 0.0

            # set the params
            self.cof = [r, s, n]
            return


        # Switch the method
        # Trust Region Reflective
        if self.fit_method == 'trf':
            bounds = (0, self.__get_fit_bounds(x, y))
            self.cof, self.cov = curve_fit(
                self._model,
                _x, _y,
                method='trf',
                sigma=self.fit_sigma,
                p0=bounds[1],
                bounds=bounds,
                **kwargs
            )

        # Levenberg-Marquardt
        elif self.fit_method == 'lm':
            self.cof, self.cov = curve_fit(
                self.model,
                _x, _y,
                method='lm',
                sigma=self.fit_sigma,
                **kwargs
            )

        else:
            raise ValueError("fit method has to be one of ['trf', 'lm']")

    def transform(self, x):
        """Transform

        Transform a given set of lag values to the theoretical variogram
        function using the actual fitting and preprocessing parameters in
        this Variogram instance

        Parameters
        ----------
        x : numpy.array
            Array of lag values to be used as model input for the fitted
            theoretical variogram model

        Returns
        -------
        numpy.array

        """
        self.preprocessing()

        # if instance is not fitted, fit it
        if self.cof is None:
            self.fit(force=True)

        # return the result
        return self.fitted_model(x)

    @property
    def fitted_model(self):
        """Fitted Model

        Returns a callable that takes a distance value and returns a
        semivariance. This model is fitted to the current Variogram
        parameters. The function will be interpreted at return time with the
        parameters hard-coded into the function code.

        Returns
        -------
        model : callable
            The current semivariance model fitted to the current Variogram
            model parameters.
        """
        if self.cof is None:
            self.fit(force=True)

        # get the pars
        cof = self.cof

        # get the function
        func = self._model

        if self._harmonize:
            code = """model = lambda x: func(x)"""
        else:
            code = """model = lambda x: func(x, %s)""" % \
               (', '.join([str(_) for _ in cof]))

        # run the code
        loc = dict(func=func)
        exec(code, loc, loc)
        model = loc['model']

        return model

    def _calc_distances(self, force=False):
        if self._dist is not None and not force:
            return

        # if self._X is of just one dimension, concat zeros.
        if self._X.ndim == 1:
            _x = np.vstack(zip(self._X, np.zeros(len(self._X))))
        else:
            _x = self._X
        # else calculate the distances
        self._dist = self._dist_func_wrapper(_x)

    def _calc_diff(self, force=False):
        """Calculates the pairwise differences

        Returns
        -------
        void

        """
        # already calculated
        if self._diff is not None and not force:
            return

        v = self.values

        # Append a column of zeros to make pdist happy
        # euclidean: sqrt((a-b)**2 + (0-0)**2) == sqrt((a-b)**2) == abs(a-b)
        self._diff = pdist(np.column_stack((v, np.zeros(len(v)))), metric="euclidean")

    def _calc_groups(self, force=False):
        """Calculate the lag class mask array

        Returns
        -------

        """
        # already calculated
        if self._groups is not None and not force:
            return

        # get the bin edges and distances
        bin_edges = self.bins
        d = self.distance

        # -1 is the group fir distances outside maxlag
        self._groups = np.ones(len(d), dtype=int) * -1

        for i, bounds in enumerate(zip([0] + list(bin_edges), bin_edges)):
            self._groups[np.where((d >= bounds[0]) & (d < bounds[1]))] = i

    def clone(self):
        """Deep copy of self

        Return a deep copy of self.

        Returns
        -------
        Variogram

        """
        return copy.deepcopy(self)

    @property
    def experimental(self):
        """Experimental Variogram

        Array of experimental (empirical) semivariance values. The array
        length will be aligned to Variogram.bins. The current
        Variogram.estimator has been used to calculate the values. Depending
        on the setting of Variogram.harmonize (True | False), either
        Variogram._experimental or Variogram.isotonic will be returned.

        Returns
        -------
        vario : numpy.ndarray
            Array of the experimental semi-variance values aligned to
            Variogram.bins.

        See Also
        --------
        Variogram._experimental
        Variogram.isotonic

        """
        return self._experimental

    @property
    def _experimental(self):
        """

        Returns
        -------

        """
        # prepare the result array
        y = np.zeros(len(self.bins), dtype=np.float64)

        # args, can set the bins for entropy
        # and should set p of percentile, not properly implemented
        if self._estimator.__name__ == 'entropy':
            bins = np.linspace(
                np.min(self.distance),
                np.max(self.distance),
                50
            )
            # apply
            for i, lag_values in enumerate(self.lag_classes()):
                y[i] = self._estimator(lag_values, bins=bins)

        # default
        else:
            for i, lag_values in enumerate(self.lag_classes()):
                y[i] = self._estimator(lag_values)

        # apply
        return y.copy()

    def __get_fit_bounds(self, x, y):
        """
        Return the bounds for parameter space in fitting a variogram model.
        The bounds are depended on the Model that is used

        Returns
        -------
        list

        """
        mname = self._model.__name__

        # use range, sill and smoothness parameter
        if mname == 'matern':
            # a is max(x), C0 is max(y) s is limited to 20?
            bounds = [np.nanmax(x), np.nanmax(y), 20.]

        # use range, sill and shape parameter
        elif mname == 'stable':
            # a is max(x), C0 is max(y) s is limited to 2?
            bounds = [np.nanmax(x), np.nanmax(y), 2.]

        # use only sill
        elif mname == 'nugget':
            # a is max(x):
            bounds = [np.nanmax(x)]

        # use range and sill
        else:
            # a is max(x), C0 is max(y)
            bounds = [np.nanmax(x), np.nanmax(y)]

        # if use_nugget is True add the nugget
        if self.use_nugget:
            bounds.append(0.99)

        return bounds

    def data(self, n=100, force=False):
        """Theoretical variogram function

        Calculate the experimental variogram and apply the binning. On
        success, the variogram model will be fitted and applied to n lag
        values. Returns the lags and the calculated semi-variance values.
        If force is True, a clean preprocessing and fitting run will be
        executed.

        Parameters
        ----------
        n : integer
            length of the lags array to be used for fitting. Defaults to 100,
            which will be fine for most plots
        force: boolean
            If True, the preprocessing and fitting will be executed as a
            clean run. This will force all intermediate results to be
            recalculated. Defaults to False

        Returns
        -------
        variogram : tuple
            first element is the created lags array
            second element are the calculated semi-variance values

        """
        # force a clean preprocessing if needed
        self.preprocessing(force=force)

        # calculate the experimental variogram
        _exp = self.experimental
        _bin = self.bins

        x = np.linspace(0, np.float64(np.nanmax(_bin)), n)

        # fit if needed
        if self.cof is None:
            self.fit(force=force)

        return x, self._model(x, *self.cof)

    @property
    def residuals(self):
        """Model residuals

        Calculate the model residuals defined as the differences between the
        experimental variogram and the theoretical model values at
        corresponding lag values

        Returns
        -------
        numpy.ndarray

        """
        # get the deviations
        experimental, model = self.model_deviations()

        # calculate the residuals
        return np.fromiter(
            map(lambda x, y: x - y, model, experimental),
            np.float
        )

    @property
    def mean_residual(self):
        """Mean Model residuals

        Calculates the mean, absoulte deviations between the experimental
        variogram and theretical model values.

        Returns
        -------
        float
        """
        return np.nanmean(np.fromiter(map(np.abs, self.residuals), np.float))

    @property
    def rmse(self):
        r"""RMSE

        Calculate the Root Mean squared error between the experimental
        variogram and the theoretical model values at corresponding lags.
        Can be used as a fitting quality measure.

        Returns
        -------
        float

        See Also
        --------
        Variogram.residuals

        Notes
        -----
        The RMSE is implemented like:

        .. math::
            RMSE = \sqrt{\frac{\sum_{i=0}^{i=N(x)} (x-y)^2}{N(x)}}

        """
        # get the deviations
        experimental, model = self.model_deviations()

        # get the sum of squares
        rsum = np.nansum(np.fromiter(
            map(lambda x, y: (x - y)**2, experimental, model),
            np.float
        ))

        return np.sqrt(rsum / len(model))

    @property
    def nrmse(self):
        r"""NRMSE

        Calculate the normalized root mean squared error between the
        experimental variogram and the theoretical model values at
        corresponding lags. Can be used as a fitting quality measure

        Returns
        -------
        float

        See Also
        --------
        Variogram.residuals
        Variogram.rmse

        Notes
        -----

        The NRMSE is implemented as:

        .. math::

            NRMSE = \frac{RMSE}{mean(y)}

        where RMSE is Variogram.rmse and y is Variogram.experimental


        """
        return self.rmse / np.nanmean(self.experimental)

    @property
    def nrmse_r(self):
        r"""NRMSE

        Alternative normalized root mean squared error between the
        experimental variogram and the theoretical model values at
        corresponding lags. Can be used as a fitting quality measure.

        Returns
        -------
        float

        See Also
        --------
        Variogram.rmse
        Variogram.nrmse

        Notes
        -----
        Unlike Variogram.nrmse, nrmse_r is not normalized to the mean of y,
        but the differece of the maximum y to its mean:

        .. math::
            NRMSE_r = \frac{RMSE}{max(y) - mean(y)}

        """
        _y = self.experimental
        return self.rmse / (np.nanmax(_y) - np.nanmean(_y))

    @property
    def r(self):
        """
        Pearson correlation of the fitted Variogram

        :return:
        """
        # get the experimental and theoretical variogram and cacluate means
        experimental, model = self.model_deviations()
        mx = np.nanmean(experimental)
        my = np.nanmean(model)

        # claculate the single pearson correlation terms
        term1 = np.nansum(np.fromiter(map(lambda x, y: (x-mx) * (y-my), experimental, model), np.float))

        t2x = np.nansum(np.fromiter(map(lambda x: (x-mx)**2, experimental), np.float))
        t2y = np.nansum(np.fromiter(map(lambda y: (y-my)**2, model), np.float))

        return term1 / (np.sqrt(t2x * t2y))

    @property
    def NS(self):
        """
        Nash Sutcliffe efficiency of the fitted Variogram

        :return:
        """
        experimental, model = self.model_deviations()
        mx = np.nanmean(experimental)

        # calculate the single nash-sutcliffe terms
        term1 = np.nansum(np.fromiter(map(lambda x, y: (x - y)**2, experimental, model), np.float))
        term2 = np.nansum(np.fromiter(map(lambda x: (x - mx)**2, experimental), np.float))

        return 1 - (term1 / term2)

    def model_deviations(self):
        """Model Deviations

        Calculate the deviations between the experimental variogram and the
        recalculated values for the same bins using the fitted theoretical
        variogram function. Can be utilized to calculate a quality measure
        for the variogram fit.

        Returns
        -------
        deviations : tuple
            first element is the experimental variogram
            second element are the corresponding values of the theoretical
            model.

        """
        # get the experimental values and their bin bounds
        _exp = self.experimental
        _bin = self.bins

        # get the model parameters
        param = self.describe()
        if 'error' in param:
            raise RuntimeError('The variogram cannot be calculated.')

        # calculate the model values at bin bounds
        _model = self.transform(_bin)

        return _exp, _model

    def describe(self):
        """Variogram parameters

        Return a dictionary of the variogram parameters.

        Returns
        -------
        dict

        """
        # fit, if not already done
        if self.cof is None:
            self.fit(force=True)

        # scale sill and range
        if self.normalized:
            maxlag = np.nanmax(self.bins)
            maxvar = np.nanmax(self.experimental)
        else:
            maxlag = 1.
            maxvar = 1.

        # get the fitting coefficents
        cof = self.cof

        # build the dict
        rdict = dict(
            name=self._model.__name__,
            estimator=self._estimator.__name__,
            effective_range=cof[0] * maxlag,
            sill=cof[1] * maxvar,
            nugget=cof[-1] * maxvar if self.use_nugget else 0
        )

        # handle s parameters for matern and stable model
        if self._model.__name__ == 'matern':
            rdict['smoothness'] = cof[2]
        elif self._model.__name__ == 'stable':
            rdict['shape'] = cof[2]

        # return
        return rdict

    @property
    def parameters(self):
        """
        Extract just the variogram parameters range, sill and nugget from the self.describe return

        :return:
        """
        d = self.describe()
        if 'error' in d:
            return [None, None, None]
        elif self._model.__name__ == 'matern':
            return list([
                d['effective_range'],
                d['sill'],
                d['smoothness'],
                d['nugget']
            ])
        elif self._model.__name__ == 'stable':
            return list([
                d['effective_range'],
                d['sill'],
                d['shape'],
                d['nugget']
            ])
        elif self._model.__name__ == 'nugget':
            return list([d['nugget']])
        else:
            return list([
                d['effective_range'],
                d['sill'],
                d['nugget']
            ])

    def to_DataFrame(self, n=100, force=False):
        """Variogram DataFrame

        Returns the fitted theoretical variogram as a pandas.DataFrame
        instance. The n and force parameter control the calaculation,
        refer to the data funciton for more info.

        Parameters
        ----------
        n : integer
            length of the lags array to be used for fitting. Defaults to 100,
            which will be fine for most plots
        force: boolean
            If True, the preprocessing and fitting will be executed as a
            clean run. This will force all intermediate results to be
            recalculated. Defaults to False

        Returns
        -------
        pandas.DataFrame

        See Also
        --------
        Variogram.data

        """
        lags, data = self.data(n=n, force=force)

        return DataFrame({
            'lags': lags,
            self._model.__name__: data}
        ).copy()

    def plot(self, axes=None, grid=True, show=True, hist=True):
        """Variogram Plot

        Plot the experimental variogram, the fitted theoretical function and
        an histogram for the lag classes. The axes attribute can be used to
        pass a list of AxesSubplots or a single instance to the plot
        function. Then these Subplots will be used. If only a single instance
        is passed, the hist attribute will be ignored as only the variogram
        will be plotted anyway.

        Parameters
        ----------
        axes : list, tuple, array, AxesSubplot or None
            If None, the plot function will create a new matplotlib figure.
            Otherwise a single instance or a list of AxesSubplots can be
            passed to be used. If a single instance is passed, the hist
            attribute will be ignored.
        grid : bool
            Defaults to True. If True a custom grid will be drawn through
            the lag class centers
        show : bool
            Defaults to True. If True, the show method of the passed or
            created matplotlib Figure will be called before returning the
            Figure. This should be set to False, when used in a Notebook,
            as a returned Figure object will be plotted anyway.
        hist : bool
            Defaults to True. If False, the creation of a histogram for the
            lag classes will be suppressed.

        Returns
        -------
        matplotlib.Figure

        """
        # get the parameters
        _bins = self.bins
        _exp = self.experimental
        x = np.linspace(0, np.nanmax(_bins), 100)  # make the 100 a param?

        # do the plotting
        if axes is None:
            if hist:
                fig = plt.figure(figsize=(8, 5))
                ax1 = plt.subplot2grid((5, 1), (1, 0), rowspan=4)
                ax2 = plt.subplot2grid((5, 1), (0, 0), sharex=ax1)
                fig.subplots_adjust(hspace=0)
            else:
                fig, ax1 = plt.subplots(1, 1, figsize=(8, 4))
                ax2 = None
        elif isinstance(axes, (list, tuple, np.ndarray)):
            ax1, ax2 = axes
            fig = ax1.get_figure()
        else:
            ax1 = axes
            ax2 = None
            fig = ax1.get_figure()

        # apply the model
        y = self.transform(x)

        # handle the relative experimental variogram
        if self.normalized:
            _bins /= np.nanmax(_bins)
            y /= np.max(_exp)
            _exp /= np.nanmax(_exp)
            x /= np.nanmax(x)

        # ------------------------
        # plot Variograms
        ax1.plot(_bins, _exp, '.b')
        ax1.plot(x, y, '-g')

        # ax limits
        if self.normalized:
            ax1.set_xlim([0, 1.05])
            ax1.set_ylim([0, 1.05])
        if grid:
            ax1.grid(False)
            ax1.vlines(_bins, *ax1.axes.get_ybound(), colors=(.85, .85, .85),
                       linestyles='dashed')
        # annotation
        ax1.axes.set_ylabel('semivariance (%s)' % self._estimator.__name__)
        ax1.axes.set_xlabel('Lag (-)')

        # ------------------------
        # plot histogram
        if ax2 is not None and hist:
            # calc the histogram
            _count = np.fromiter(
                (g.size for g in self.lag_classes()), dtype=int
            )

            # set the sum of hist bar widths to 70% of the x-axis space
            w = (np.max(_bins) * 0.7) / len(_count)

            # plot
            ax2.bar(_bins, _count, width=w, align='center', color='red')

            # adjust
            plt.setp(ax2.axes.get_xticklabels(), visible=False)
            ax2.axes.set_yticks(ax2.axes.get_yticks()[1:])

            # need a grid?
            if grid:  #pragma: no cover
                ax2.grid(False)
                ax2.vlines(_bins, *ax2.axes.get_ybound(),
                           colors=(.85, .85, .85), linestyles='dashed')

            # anotate
            ax2.axes.set_ylabel('N')

        # show the figure
        if show:  # pragma: no cover
            fig.show()

        return fig

    def scattergram(self, ax=None, show=True):

        # create a new plot or use the given
        if ax is None:
            fig, ax = plt.subplots(1, 1)
        else:
            fig = ax.get_figure()

        tail = np.empty(0)
        head = tail.copy()

        for h in np.unique(self.lag_groups()):
            # get the head and tail
            x, y = np.where(squareform(self.lag_groups()) == h)

            # concatenate
            tail = np.concatenate((tail, self.values[x]))
            head = np.concatenate((head, self.values[y]))

        # plot the mean on tail and head
        ax.vlines(np.mean(tail), np.min(tail), np.max(tail), linestyles='--',
                  color='red', lw=2)
        ax.hlines(np.mean(head), np.min(head), np.max(head), linestyles='--',
                  color='red', lw=2)
        # plot
        ax.scatter(tail, head, 10, marker='o', color='orange')

        # annotate
        ax.set_ylabel('head')
        ax.set_xlabel('tail')

        # show the figure
        if show:  # pragma: no cover
            fig.show()

        return fig

    def location_trend(self, axes=None, show=True):
        """Location Trend plot

        Plots the values over each dimension of the coordinates in a scatter
        plot. This will visually show correlations between the values and any
        of the coordinate dimension. If there is a value dependence on the
        location, this would violate the intrinsic hypothesis. This is a
        weaker form of stationarity of second order.

        Parameters
        ----------
        axes : list
            Can be None (default) or a list of matplotlib.AxesSubplots. If a
            list is passed, the location trend plots will be plotted on the
            given instances. Note that then length of the list has to match
            the dimeonsionality of the coordinates array. In case 3D
            coordinates are used, three subplots have to be given.

        Returns
        -------
        matplotlib.Figure

        """
        N = len(self._X[0])
        if axes is None:
            # derive the needed amount of col and row
            nrow = int(round(np.sqrt(N)))
            ncol = int(np.ceil(N / nrow))
            fig, axes = plt.subplots(nrow, ncol, figsize=(ncol * 6 ,nrow * 6))
        else:
            if not len(axes) == N:
                raise ValueError(
                    'The amount of passed axes does not fit the coordinate' +
                    ' dimensionality of %d' % N)
            fig = axes[0].get_figure()

        for i in range(N):
            axes.flatten()[i].plot([_[i] for _ in self._X], self.values, '.r')
            axes.flatten()[i].set_xlabel('%d-dimension' % (i + 1))
            axes.flatten()[i].set_ylabel('value')

        # plot the figure and return it
        plt.tight_layout()

        if show:  # pragma: no cover
            fig.show()

        return fig

    def distance_difference_plot(self, ax=None, plot_bins=True, show=True):
        """Raw distance plot

        Plots all absoulte value differences of all point pair combinations
        over their separating distance, without sorting them into a lag.

        Parameters
        ----------
        ax : None, AxesSubplot
            If None, a new matplotlib.Figure will be created. In case a
            Figure was already created, pass the Subplot to use as ax argument.
        plot_bins : bool
            If True (default) the bin edges will be included into the plot.
        show : bool
            If True (default), the show method of the Figure will be called
            before returning the Figure. Can be set to False, to avoid
            doubled figure rendering in Jupyter notebooks.

        Returns
        -------
        matplotlib.pyplot.Figure

        """
        # get all distances
        _dist = self.distance

        # get all differences
        if self._diff is None:
            self._calc_diff()
        _diff = self._diff

        # create the plot
        if ax is None:
            fig, ax = plt.subplots(1, 1, figsize=(8, 6))
        else:
            fig = ax.get_figure()

        # plot the bins
        if plot_bins:
            _bins = self.bins
            ax.vlines(_bins, 0, np.max(_diff), linestyle='--', lw=1, color='r')

        # plot
        ax.scatter(_dist, _diff, 8, color='b', marker='o', alpha=0.5)

        # set limits
        ax.set_ylim((0, np.max(_diff)))
        ax.set_xlim((0, np.max(_dist)))
        ax.set_xlabel('separating distance')
        ax.set_ylabel('pairwise difference')
        ax.set_title('Pairwise distance ~ difference')

        # show the plot
        if show:  # pragma: no cover
            fig.show()

        return fig

    def __repr__(self):  # pragma: no cover
        """
        Textual representation of this Variogram instance.

        :return:
        """
        try:
            _name = self._model.__name__
            _b = int(len(self.bins))
        except:
            return "< abstract Variogram >"
        return "< %s Semivariogram fitted to %d bins >" % (_name, _b)

    def __str__(self):  # pragma: no cover
        """String Representation

        Descriptive respresentation of this Variogram instance that shall give
        the main variogram parameters in a print statement.

        Returns
        -------
        description : str
            String description of the variogram instance. Described by the
            Variogram parameters.

        """
        par = self.describe()

        _sill = np.NaN if 'error' in par else par['sill']
        _range = np.NaN if 'error' in par else par['effective_range']
        _nugget = np.NaN if 'error' in par else par['nugget']

        s = "{0} Variogram\n".format(par['name'])
        s+= "-" * (len(s) - 1) + "\n"
        s+="""Estimator:         %s
        \rEffective Range:   %.2f
        \rSill:              %.2f
        \rNugget:            %.2f
        """ % (par['estimator'], _range, _sill, _nugget)

        return s