docs: def _initialize_variogram_model(X, y, variogram_model,
                                  variogram_model_parameters, variogram_function,
                                  nlags, weight, coordinates_type):
      """Initializes the variogram model for kriging. If user does not specify
      parameters, calls automatic variogram estimation routine.
      Returns lags, semivariance, and variogram model parameters.
  
      Parameters
      ----------
      X: ndarray
          float array [n_samples, n_dim], the input array of coordinates
      y: ndarray
          float array [n_samples], the input array of values to be kriged
      variogram_model: str
          user-specified variogram model to use
      variogram_model_parameters: list
          user-specified parameters for variogram model
      variogram_function: callable
          function that will be called to evaluate variogram model
          (only used if user does not specify variogram model parameters)
      nlags: int
          integer scalar, number of bins into which to group inter-point distances
      weight: bool
          boolean flag that indicates whether the semivariances at smaller lags
          should be weighted more heavily in the automatic variogram estimation
      coordinates_type: str
          type of coordinates in X array, can be 'euclidean' for standard
          rectangular coordinates or 'geographic' if the coordinates are lat/lon
  
      Returns
      -------
      lags: ndarray
          float array [nlags], distance values for bins into which the
          semivariances were grouped
      semivariance: ndarray
          float array [nlags], averaged semivariance for each bin
      variogram_model_parameters: list
          parameters for the variogram model, either returned unaffected if the
          user specified them or returned from the automatic variogram
          estimation routine
      """
  
      # distance calculation for rectangular coords now leverages
      # scipy.spatial.distance's pdist function, which gives pairwise distances
      # in a condensed distance vector (distance matrix flattened to a vector)
      # to calculate semivariances...
      if coordinates_type == 'euclidean':
          d = pdist(X, metric='euclidean')
          g = 0.5 * pdist(y[:, None], metric='sqeuclidean')
  
      # geographic coordinates only accepted if the problem is 2D
      # assume X[:, 0] ('x') => lon, X[:, 1] ('y') => lat
      # old method of distance calculation is retained here...
      # could be improved in the future
      elif coordinates_type == 'geographic':
          if X.shape[1] != 2:
              raise ValueError('Geographic coordinate type only '
                               'supported for 2D datasets.')
          x1, x2 = np.meshgrid(X[:, 0], X[:, 0], sparse=True)
          y1, y2 = np.meshgrid(X[:, 1], X[:, 1], sparse=True)
          z1, z2 = np.meshgrid(y, y, sparse=True)
          d = great_circle_distance(x1, y1, x2, y2)
          g = 0.5 * (z1 - z2)**2.
          indices = np.indices(d.shape)
          d = d[(indices[0, :, :] > indices[1, :, :])]
          g = g[(indices[0, :, :] > indices[1, :, :])]
  
      else:
          raise ValueError("Specified coordinate type '%s' "
                           "is not supported." % coordinates_type)
  
      # Equal-sized bins are now implemented. The upper limit on the bins
      # is appended to the list (instead of calculated as part of the
      # list comprehension) to avoid any numerical oddities
      # (specifically, say, ending up as 0.99999999999999 instead of 1.0).
      # Appending dmax + 0.001 ensures that the largest distance value
      # is included in the semivariogram calculation.
      dmax = np.amax(d)
      dmin = np.amin(d)
      dd = (dmax - dmin) / nlags
      bins = [dmin + n * dd for n in range(nlags)]
      dmax += 0.001
      bins.append(dmax)
  
      # This old binning method was experimental and doesn't seem
      # to work too well. Bins were computed such that there are more
      # at shorter lags. This effectively weights smaller distances more
      # highly in determining the variogram. As Kitanidis points out,
      # the variogram fit to the data at smaller lag distances is more
      # important. However, the value at the largest lag probably ends up
      # being biased too high for the larger values and thereby throws off
      # automatic variogram calculation and confuses comparison of the
      # semivariogram with the variogram model.
      #
      # dmax = np.amax(d)
      # dmin = np.amin(d)
      # dd = dmax - dmin
      # bins = [dd*(0.5**n) + dmin for n in range(nlags, 1, -1)]
      # bins.insert(0, dmin)
      # bins.append(dmax)
  
      lags = np.zeros(nlags)
      semivariance = np.zeros(nlags)
  
      for n in range(nlags):
          # This 'if... else...' statement ensures that there are data
          # in the bin so that numpy can actually find the mean. If we
          # don't test this first, then Python kicks out an annoying warning
          # message when there is an empty bin and we try to calculate the mean.
          if d[(d >= bins[n]) & (d < bins[n + 1])].size > 0:
              lags[n] = np.mean(d[(d >= bins[n]) & (d < bins[n + 1])])
              semivariance[n] = np.mean(g[(d >= bins[n]) & (d < bins[n + 1])])
          else:
              lags[n] = np.nan
              semivariance[n] = np.nan
  
      lags = lags[~np.isnan(semivariance)]
      semivariance = semivariance[~np.isnan(semivariance)]
  
      # a few tests the make sure that, if the variogram_model_parameters
      # are supplied, they have been supplied as expected...
      # if variogram_model_parameters was not defined, then estimate the variogram
      if variogram_model_parameters is not None:
          if variogram_model == 'linear' and len(variogram_model_parameters) != 2:
              raise ValueError("Exactly two parameters required "
                               "for linear variogram model.")
          elif variogram_model in ['power', 'spherical', 'exponential',
                                   'gaussian', 'hole-effect'] \
                  and len(variogram_model_parameters) != 3:
              raise ValueError("Exactly three parameters required for "
                               "%s variogram model" % variogram_model)
      else:
          if variogram_model == 'custom':
              raise ValueError("Variogram parameters must be specified when "
                               "implementing custom variogram model.")
          else:
              variogram_model_parameters = \
                  _calculate_variogram_model(lags, semivariance, variogram_model,
                                             variogram_function, weight)
  
      return lags, semivariance, variogram_model_parameters