docs: def _calculate_variogram_model(lags, semivariance, variogram_model,
                                 variogram_function, weight):
      """Function that fits a variogram model when parameters are not specified.
      Returns variogram model parameters that minimize the RMSE between the
      specified variogram function and the actual calculated variogram points.
  
      Parameters
      ----------
      lags: 1d array
          binned lags/distances to use for variogram model parameter estimation
      semivariance: 1d array
          binned/averaged experimental semivariances to use for variogram model
          parameter estimation
      variogram_model: str/unicode
          specified variogram model to use for parameter estimation
      variogram_function: callable
          the actual funtion that evaluates the model variogram
      weight: bool
          flag for implementing the crude weighting routine, used in order to fit
          smaller lags better this is passed on to the residual calculation
          cfunction, where weighting is actually applied...
  
      Returns
      -------
      res: list
          list of estimated variogram model parameters
  
      NOTE that the estimation routine works in terms of the partial sill
      (psill = sill - nugget) -- setting bounds such that psill > 0 ensures that
      the sill will always be greater than the nugget...
      """
  
      if variogram_model == 'linear':
          x0 = [(np.amax(semivariance) - np.amin(semivariance)) /
                (np.amax(lags) - np.amin(lags)), np.amin(semivariance)]
          bnds = ([0., 0.], [np.inf, np.amax(semivariance)])
      elif variogram_model == 'power':
          x0 = [(np.amax(semivariance) - np.amin(semivariance)) /
                (np.amax(lags) - np.amin(lags)), 1.1, np.amin(semivariance)]
          bnds = ([0., 0.001, 0.], [np.inf, 1.999, np.amax(semivariance)])
      else:
          x0 = [np.amax(semivariance) - np.amin(semivariance),
                0.25*np.amax(lags), np.amin(semivariance)]
          bnds = ([0., 0., 0.], [10.*np.amax(semivariance), np.amax(lags),
                                 np.amax(semivariance)])
  
      # use 'soft' L1-norm minimization in order to buffer against
      # potential outliers (weird/skewed points)
      res = least_squares(_variogram_residuals, x0, bounds=bnds, loss='soft_l1',
                          args=(lags, semivariance, variogram_function, weight))
  
      return res.x