docs: def minimizeCompass(func, x0, args=(),
              bounds=None, scaling=None,
              redfactor=2.0, deltainit=1.0, deltatol=1e-3, feps=1e-15,
              errorcontrol=True, funcNinit=30, funcmultfactor=2.0,
              paired=True, alpha=0.05, disp=False, callback=None, **kwargs):
      """
      Minimization of an objective function by a pattern search.
  
      The algorithm does a compass search along coordinate directions.
      If `errorcontrol=True` then the function is called repeatedly to average
      over the stochasticity in the function evaluation. The number of
      evaluations over which to average is adapted dynamically to ensure convergence.
  
      The algorithm terminates when the current iterate is locally optimally
      at the target pattern size deltatol or when the function value differs by
      less than the tolerance feps along all directions. 
  
      Parameters
      ----------
      func: callable
          objective function to be minimized:
          called as `func(x, *args)`,
          if `paired=True`, then called with keyword argument `seed` additionally
      x0: array-like
          starting point
      args: tuple
          extra arguments to be supplied to func
      bounds: array-like
          bounds on the variables
      scaling: array-like
          scaling by which to multiply step size and tolerances along different dimensions
      redfactor: float
          reduction factor by which to reduce delta if no reduction direction found 
      deltainit: float
          initial pattern size
      deltatol: float
          target pattern size,
          function differences at this scale need to be larger than stochasticitiy
          in evaluations to ensure convergence if `errorcontrol=False`
      feps: float
         smallest difference in function value to resolve 
      errorcontrol: boolean
          whether to control error of simulation by repeated sampling
      funcNinit: int, only for errorcontrol=True
          initial number of iterations to use for the function, do not set much lower
          than 30 as otherwise there is no sufficient statistics for function comparisons
      funcmultfactor: float, only for errorcontrol=True
          multiplication factor by which to increase number of iterations of function
      paired: boolean, only for errorcontrol=True
          compare for same random seeds
      alpha: float, only for errorcontrol=True
          significance level of tests, the higher this value the more statistics
          is acquired, which decreases the risk of taking a step in a non-descent
          direction at the expense of higher computational cost per iteration
      disp: boolean
          whether to output status updates during the optimization
      callback: callable
          called after each iteration, as callback(xk), where xk is the current parameter vector.
  
      Returns
      -------
      scipy.optimize.OptimizeResult object
          special entry: free
          Boolean array indicating parameters that are unconstrained at the optimum (within feps) 
      """
      #TODO: implement variable deltas for different directions (might speed up things, see review)
      if disp:
          print('minimization starting')
          print('args', args)
          print('errorcontrol', errorcontrol)
          print('paired', paired)
      # absolute tolerance for float comparisons
      floatcompatol = 1e-14
      x0 = np.asarray(x0)
      if scaling is None:
          scaling = np.ones(x0.shape)
      else:
          scaling = np.asarray(scaling)
      # ensure initial point lies within bounds
      if bounds is not None:
          bounds = np.asarray(bounds)
          np.clip(x0, bounds[:, 0], bounds[:, 1], out=x0)
      def clip(x, d):
          """clip x+d to respect bounds
          returns clipped result and effective distance"""
          xnew = x + d
          if bounds is not None:
              # if test point depasses set to boundary instead
              xclipped = np.clip(xnew, bounds[:, 0], bounds[:, 1])
              deltaeff = np.abs(x - xclipped).sum()
              return xclipped, deltaeff
          else:
              return xnew, delta
      # generate set of search directions (+- s_i e_i | i = 1, ...,  N)
      def unit(i, N):
          "return ith unit vector in R^N"
          arr = np.zeros(N)
          arr[i] = 1.0
          return arr
      N = len(x0)
      generatingset = [unit(i, N)*direction*scaling[i] for i in np.arange(N) for direction in [+1, -1]]
     
      # memoize function
      if errorcontrol:
          funcm = AveragedFunction(
              func, fargs=args, paired=paired, N=funcNinit)
          # apply Bonferroni correction to confidence level
          # (need more statistics in higher dimensions)
          alpha = alpha/len(generatingset)
      else:
          # freeze function arguments
          def funcf(x, **kwargs):
              return func(x, *args, **kwargs)
          funcm = _memoized(funcf)
  
      x = x0 
      delta = deltainit
      # number of iterations
      nit = 0
      # continue as long as delta is larger than tolerance
      # or if there was an update during the last iteration
      found = False
      while delta >= deltatol-floatcompatol or found:
          nit += 1
          # if delta gets close to deltatol, do iteration with step size deltatol instead
          if delta/redfactor < deltatol:
              delta = deltatol
          if disp:
              print('nit %i, Delta %g' % (nit, delta))
          found = False
          np.random.shuffle(generatingset)
          for d in generatingset:
              xtest, deltaeff = clip(x, delta*d)
              if deltaeff < floatcompatol:
                  continue
              # Does xtest improve upon previous function value?
              if ((not errorcontrol and (funcm(xtest) < funcm(x)-feps))
                 or (errorcontrol
                     and funcm.test(xtest, x, type_='smaller', alpha=alpha))):
                  x = xtest
                  found = True
                  if disp:
                      print(x)
              # Is non-improvement due to too large step size or missing statistics?
              elif ((deltaeff >= deltatol*np.sum(np.abs(d))) # no refinement for boundary steps smaller than tolerance
                      and ((not errorcontrol and (funcm(xtest) < funcm(x)+feps))
                          or (errorcontrol
                              and funcm.test(xtest, x, type_='equality', alpha=alpha)
                              and (funcm.diffse(xtest, x) > feps)))):
                  # If there is no significant difference the step size might
                  # correspond to taking a step to the other side of the minimum.
                  # Therefore test if middle point is better
                  xmid = 0.5*(x+xtest)
                  if ((not errorcontrol and funcm(xmid) < funcm(x)-feps)
                      or (errorcontrol
                          and funcm.test(xmid, x, type_='smaller', alpha=alpha))):
                      x = xmid
                      delta /= redfactor
                      found = True
                      if disp:
                          print('mid', x)
                  # otherwise increase accuracy of simulation to try to get to significance
                  elif errorcontrol:
                      funcm.N *= funcmultfactor
                      if disp:
                          print('new N %i' % funcm.N)
                      found = True
          if callback is not None:
              callback(x)
          if not found:
              delta /= redfactor
  
      message = 'convergence within deltatol'
      # check if any of the directions are free at the optimum
      delta = deltatol
      free = np.zeros(x.shape, dtype=bool)
      for d in generatingset:
          dim = np.argmax(np.abs(d))
          xtest, deltaeff = clip(x, delta*d)
          if deltaeff < deltatol*np.sum(np.abs(d))-floatcompatol: # do not consider as free for boundary steps
              continue
          if not free[dim] and (((not errorcontrol and funcm(xtest) - feps < funcm(x)) or
              (errorcontrol and funcm.test(xtest, x, type_='equality', alpha=alpha)
                  and (funcm.diffse(xtest, x) < feps)))):
              free[dim] = True
              message += '. dim %i is free at optimum' % dim
                  
      reskwargs = dict(x=x, nit=nit, nfev=funcm.nev, message=message, free=free,
                       success=True)
      if errorcontrol:
          f, funse = funcm(x)
          res = OptimizeResult(fun=f, funse=funse, **reskwargs)
      else:
          f = funcm(x)
          res = OptimizeResult(fun=f, **reskwargs)
      if disp:
          print(res)
      return res