docs: def clean_retinotopy(hemi, retinotopy=Ellipsis, mask=Ellipsis, weight=Ellipsis,
                       surface='midgray', min_weight=Ellipsis, min_eccentricity=0.75,
                       visual_area=Ellipsis, map_visual_areas=Ellipsis,
                       visual_area_field_signs=Ellipsis,
                       measurement_uncertainty=0.3, measurement_knob=1,
                       magnification_knob=0, fieldsign_knob=6, edge_knob=0,
                       yield_report=False, steps=100, rounds=4, output_style='visual',
                       jitter=None, average=None):
  
      '''
      clean_retinotopy(hemi) attempts to cleanup the retinotopic maps on the given cortical mesh by
        minimizing an objective function that tracks the smoothness of the fields, the orthogonality
        of polar angle to eccentricity, and the deviation of the values from the measured values; the
        yielded result is the smoothed retinotopy, as would be returned by
        as_retinotopy(..., 'visual').
  
      The argument hemi may be a Cortex object or a Mesh. See also clean_retinotopy_potential for
      information on the various parameters related to the potential function that is minimized by
      clean_retinotopy(). The following additional options are also accepted:
        * output_style (default: 'visual') specifies the style of the output data that should be
          returned; this should be a string understood by as_retinotopy.
        * yield_report (False) may be set to True, in which case a tuple (retino, report) is returned,
          where the report is the return value of the scipy.optimization.minimize function.
      '''
      # parse our args
      if jitter in [True, Ellipsis, 'auto', 'automatic']: jitter = (4, 0.05, 1)
      if is_tuple(jitter) and len(jitter) > 0:
          if len(jitter) > 3:  raise ValueError('jitter tuple must be (mod, scale, phase)')
          if len(jitter) == 1: jitter = jitter + (0.005,)
          if len(jitter) == 2: jitter = jitter + (1,)
          (jitter_mod, jitter_scale, jitter_phase) = jitter
      else: jitter = None
      if average in [True, Ellipsis, 'auto', 'automatic']: average = (4, 3)
      if is_tuple(average) and len(average) > 0:
          if len(average) > 2:  raise ValueError('average tuple must be (mod, phase)')
          if len(average) == 1: average = average + (3,)
          (average_mod, average_phase) = average
      else: average = None
      if visual_area_field_signs is None: visual_area_field_signs = {}
      # First, make the potential function:
      f = clean_retinotopy_potential(hemi, retinotopy=retinotopy, mask=mask, weight=weight,
                                     surface=surface, min_weight=min_weight,
                                     min_eccentricity=min_eccentricity,
                                     measurement_uncertainty=measurement_uncertainty,
                                     measurement_knob=measurement_knob,
                                     magnification_knob=magnification_knob,
                                     fieldsign_knob=fieldsign_knob, edge_knob=edge_knob,
                                     visual_area=visual_area, map_visual_areas=map_visual_areas,
                                     visual_area_field_signs=visual_area_field_signs)
      # at this point, it's possible that we got a lazy map back; if so we're going to want to iterate
      # through it; otherwise, we'll want to just iterate through the single return value...
      m = f if pimms.is_map(f) else {None: f}
      (x,y) = np.full((2, hemi.vertex_count), np.nan) # the output x/y prf centers
      tess = hemi if geo.is_tess(hemi) else hemi.tess
      for (k,f) in six.iteritems(m):
          # The initial parameter vector is stored in the meta-data:
          X0 = f.meta_data['X0']
          submesh = f.meta_data['mesh']
          X = X0
          for ii in range(rounds):
              mtd = 'L-BFGS-B' if (ii % 2) == 0 else 'TNC'
              if jitter is not None and ii % jitter_mod == jitter_phase:
                  ec = np.sqrt(np.sum(X**2, axis=1))
                  th = (np.random.rand(len(ec)) - 0.5)*2*np.pi
                  r  = np.random.exponential(ec*jitter_scale)
                  X = X + np.transpose([r*np.cos(th), r*np.sin(th)])
              if average is not None and ii % average_mod == average_phase:
                  X = np.array([X[k] if len(nn) == 0 else np.mean(X[list(nn)],0)
                                for (k,nn) in enumerate(submesh.tess.indexed_neighborhoods)])
              rr = f.minimize(X, method=mtd, options=dict(maxiter=steps, disp=False))
              X = rr.x
          X = np.reshape(X, X0.shape)
          if X.shape[0] != 2: X = X.T
          for (u,v) in zip([x,y], X): u[tess.index(submesh.labels)] = v
      return as_retinotopy({'x':x, 'y':y}, output_style)