docs: def clean_retinotopy_potential(hemi, retinotopy=Ellipsis, mask=Ellipsis, weight=Ellipsis,
                                 surface='midgray', min_weight=Ellipsis, min_eccentricity=0.75,
                                 visual_area=None, 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):
      '''
      clean_retinotopy_potential(hemi) yields a retinotopic potential function for the given
        hemisphere that, when minimized, should yeild a cleaned/smoothed version of the retinotopic
        maps.
  
      The potential function f returned by clean_retinotopy_potential() is a PotentialFunction object,
      as defined in neuropythy.optimize. The potential function consists of four terms that are summed
      with weights derived from the four '*_knob' options (see below). The function f as well as the
      three terms that it comprises require as input a matrix X of the pRF centers of mesh or the
      masked part of the mesh (X0 is the initial measurement matrix). These four potential terms are:
        * The measurement potential. The retinotopic map that is being minimized is referred to as the
          measured map, and the measurement potential function, fm(X), increases as X becomes farther
          from X0. Explicitly, fm(X) is the sum over all pRF centers (x,y) in X (with initial position
          (x0,y0) in X0) of exp(-0.5 * ((x - x0)**2 + (y - y0)**2) / s**2). The parameter s is the
          initial eccentricity (sqrt(x0**2 + y0**2)) times the measurement_uncertainty option.
        * The magnification potential. The retinotopy cleaning is performed in part by attempting to
          smooth the visual magnification (the inverse of cortical magnification: deg**2 / mm**2)
          across the cortical surface; the magnification potential fg(X) specifies how the visual
          magnification contributes to the overall potential: it decreases as the magnification
          becomes smoother and increases as it becomes less smooth. Explicitly, fg(X) is equal to the
          sum over all pairs of faces (s,t) sharing one edge e of
          (vmag(s) - sgn(vmag(t))*vmag(e))**2 + (vmag(t) - sgn(vmag(s))*vmag(e))**2. Note that the
          function vmag(s) yields the areal visual magnification (deg**2 / mm**2) of any face s and
          vmag(e) is the square of the linear magnification of any edge e; additionally, the sign of
          vmag(s) for a face s is always equal to the fieldsign of the face (while for edges vmag(e)
          is always positive).
        * The fieldsign potential. The next way in which the potential function attempts to clean the
          retinotopy is via the use of fieldsign: adjacent faces should have the same fieldsign under
          most circumstanced. This is modeled by the function fs(X), which is 0 for any pair of faces
          that have the same fieldsign and non-zero for faces that have different fieldsigns. The
          form of fs(X) is the sum over all pairs of adjacent triangles (s,t) of -vmag(s)*vmag(t) if
          vmag(s) and vmag(t) have different fieldsigns, otherwise 0.
        * The edge potential. Finally, the potential function attempts to force edges to be smooth by
          penalizing edges whose endpoints are far apart in the visual field. The edge potential
          function fe(X) is equal to the sum for all edges (u,v) of
          (x[u] - x[v])**2 + (y[u] - y[v])**2 / mean(eccen(u), eccen(v)).
  
      Note additionally that all four potential functions are normalized by a factor intended to keep
      them on similar scales (this factor is not mentioned above or below, but it is automatically
      applied to all potential terms). For the magnification, fieldsign, and edge potential terms, the
      normalization factor is 1/m where m is the number of non-perimeter edges (or, alternately, the
      number of adjacent face pairs) in the mesh. For the measurement potential term, the
      normalization factor is 1/W where W is the sum of the weights on the measurement vertices (if
      no weights are given, they are considered to be 1 for each vertex).
  
      The following options may be given:
        * retinotopy (default: Ellipsis) specifies the retinotopy data to use for the hemisphere;
          the argument may be a map from retinotopy_data or a valid argument to it. The default
          indicates that the result of calling retinotopy_data(hemi) is used.
        * mask (default: Ellipsis) specifies that the specific mask should be used; by default, the
          mask is made using the vertices kept in to_flatmap('occipital_pole', hemi, radius=pi/2.75).
        * weight (default: Ellipsis) specifies the weight to use; the default indicates that the
          weight found in the retinotopy data, if any, should be used. If None, then all values
          in the mask are equally weighted.
        * visual_area (default: Ellipsis) specifies the visual area labels to use; the default
          indicates that the visual area property found in the retinotopy data should be used, if any.
          If None then no visual area splitting is done. This property is only important if 
          map_visual_areas is not False or None; otherwise it is ignored.
        * map_visual_areas (default: Ellipsis) specifies whether the return value should be a lazy map
          whose keys are visual area labels and whose values are recursed calls to this function for
          only the subset of the mesh with the associated label. May be False or None to specify that
          a single potential should be yielded. May be a list of labels to specify that only those
          visual areas should be mapped; the default value (Ellipsis) uses all labels in visual_areas
          except for 0.
        * min_weight (default: Ellipsis) specifies the minimum weight to include, after the
          weights have been normalized such that sum(weights) == 1. If the value is a list or
          tuple containing a single item [p] then p is taken to be a percentile below which
          vertices should be excluded. The default, Ellipsis, is equivalent to [5].
        * min_eccentricity (default: 0.75) specifies the eccentricity below which no measurement-based
          potential is applied; i.e., by default, vertices with eccentricity below 0.75 degrees will
          be considered as having 0 weight.
        * surface (default: 'midgray') specifies which surface should be used to establish cortical
          magnification; may be 'pial', 'midgray', or 'white'.
        * measurement_uncertainty (default: 0.3) is used to determine the standard deviation of the 
          Gaussian potential well used to prevent individual vertices with valid retinotopic
          measurements from straying too far from their initial measured positions. In other words, if
          a vertex has a weight that is above threshold and a pRF center of (x0,y0), then the 
          measurement-potential for that vertex is exp(-0.5 * ((x - x0)**2 + (y - y0)**2)/s**2) where
          (x,y) is the center of the pRF during minimization and s is equal to
          measurement_uncertainty * sqrt(x0**2 + y0**2).
        * measurement_knob, magnification_knob, fieldsign_knob, and edge_knob (defaults: 1, 0, 12, 0,
          respectively) specify the relative weights of the terms of the potential function on a log2
          scale. In other words, if the measurement, magnification, fieldsign, and edge potential
          terms are fm, fg, fs, and fe while the knobs are km, kg, ks, and ke, then the overall
          potential function f is equal to:
          f(X) = (2**km * fm(X) + 2**kg * fg(X) + 2**ks * fs(X) + 2**ke * fe(X)) / q
          where w = (2**km + 2**kg + 2**ks + 2**ke)
          If any knob is set to None, then its value is 0 instead of 2**k.
      '''
      from neuropythy.util import curry
      import neuropythy.optimize as op
      # first, get the mesh and the retinotopy data
      mesh = geo.to_mesh((hemi, surface))
      rdat = (retinotopy_data(mesh) if retinotopy is Ellipsis   else
              retinotopy            if pimms.is_map(retinotopy) else
              retinotopy_data(mesh, retinotopy))
      lbls = (rdat.get('visual_area') if visual_area is Ellipsis else
              None                    if visual_area is None     else
              hemi.property(visual_area))
      wght = (rdat.get('variance_explained') if weight is Ellipsis                      else
              weight                         if pimms.is_vector(weight)                 else
              rdat.get(weight)               if pimms.is_str(weight) and weight in rdat else
              hemi.property(weight)          if weight is not None                      else
              None)
      # figure out the mask
      if mask is Ellipsis:
          if mesh.coordinates.shape[0] == 2: mask = mesh.indices
          else: mask = geo.to_flatmap('occipital_pole', hemi, radius=np.pi/2.75).labels
      else: mask = hemi.mask(mask, indices=True)
      global_field_sign = None
      # if we are splitting on visual area, we should do that here:
      if map_visual_areas and lbls is not None:
          # get the visual areas
          vas = (np.unique(lbls)                    if map_visual_areas == 'all'           else
                 np.setdiff1d(np.unique(lbls), [0]) if map_visual_areas in [True,Ellipsis] else
                 np.unique(map_visual_areas))
          # we also want to have the field-sign map handy if provided
          if visual_area_field_signs is None: visual_area_field_signs = {}
          elif visual_area_field_signs is Ellipsis: visual_area_field_signs = {1:-1, 2:1, 3:-1, 4:1}
          # special case when map_visual_areas is an integer/string (label)
          if pimms.is_int(map_visual_areas) or pimms.is_str(map_visual_areas):
              mask = np.intersect1d(mask, np.where(lbls == map_visual_areas)[0])
              global_field_sign = visual_area_field_signs.get(map_visual_areas)
          else: # otherwise we return a lazy map
              kw = dict(retinotopy=rdat, mask=mask, weight=wght,
                        surface=surface, min_weight=min_weight, min_eccentricity=min_eccentricity,
                        visual_area=lbls, measurement_uncertainty=measurement_uncertainty,
                        measurement_knob=measurement_knob,
                        magnification_knob=magnification_knob, fieldsign_knob=fieldsign_knob,
                        edge_knob=edge_knob, visual_area_field_signs=visual_area_field_signs)
              return pimms.lazy_map(
                  {lbl: curry(clean_retinotopy_potential, hemi, map_visual_areas=lbl, **kw)
                   for lbl in vas})
      # fix rdat, weight, and mesh to match the mask
      (supermesh, orig_mask) = (mesh, mask)
      rdat = {k:(v[mask] if len(v) > len(mask) else v) for (k,v) in six.iteritems(rdat)}
      mesh = supermesh.submesh(mask)
      # possible that the mask got further downsampled:
      mask = supermesh.tess.index(mesh.labels)
      if len(mask) == len(orig_mask): smsk = np.arange(len(mask))
      else:
          tmp = set(mask)
          smsk = np.asarray([k for (k,u) in enumerate(orig_mask) if u in tmp])
      n = mesh.vertex_count # number vertices
      N = 2*n # number parameters
      if wght is None:                  wght = np.ones(n)
      elif len(wght) == len(orig_mask): wght = np.array(wght)[smsk]
      elif len(wght) > n:               wght = np.array(wght)[mask]
      else:                             wght = np.array(wght)
      wght[~np.isfinite(wght)] = 0
      if min_eccentricity < 0 or np.isclose(min_eccentricity, 0):
          raise ValueError('min_eccentricity should be small but must be > 0')
      # we'll need a potential function...
      # The input to the potential function will be a 2 x N matrix of (x,y) visual field coordinates:
      xy      = op.identity
      (x,y)   = (xy[np.arange(0,N,2)],xy[np.arange(1,N,2)])
      # We have a few components to the potential function
      # [1] Deviation from measurements:
      # These are the initial measurements we will use
      xy0     = np.array(as_retinotopy(rdat, 'geographical')).T
      if len(xy0) == len(orig_mask): xy0 = xy0[smsk]
      elif len(xy0) > n:             xy0 = xy0[mask]
      (x0,y0) = xy0.T
      xy0     = xy0.flatten()
      ecc0    = np.sqrt(x0**2 + y0**2)
      ii      = np.where(ecc0 > min_eccentricity)[0]
      minw    = (0                             if min_weight is None          else
                 np.percentile(wght[ii], 5)    if min_weight is Ellipsis      else
                 min_weight                    if pimms.is_number(min_weight) else
                 0                             if np.std(wght[ii]) < 0.00001  else
                 np.percentile(wght[ii], min_weight[0]))
      ii      = np.intersect1d(ii, np.where(wght > minw)[0])
      wsum    = np.sum(wght[ii])
      if wsum < 0 or np.isclose(wsum, 0): raise ValueError('all-zero weights given')
      wght    = wght / wsum
      s2_meas = (measurement_uncertainty * ecc0[ii])**2
      d2_meas = (x[ii] - x0[ii])**2 + (y[ii] - y0[ii])**2
      f_meas  = (1 - op.exp(-0.5*d2_meas/s2_meas)) * wght[ii]
      f_meas  = op.sum(f_meas)
      # [2] For adjacent triangles, how different are the cortical magnifications?
      sarea   = mesh.face_areas
      faces   = mesh.tess.indexed_faces.T
      selen   = mesh.edge_lengths
      # we want to ensure that vmag is locally smooth across triangles, but we need
      # to make sure there aren't any edges or faces with 0 surface-areas (so that
      # we don't divide by zero)
      mnden   = 0.0001
      (e,s,t) = np.transpose([(i,e[0],e[1]) for (i,e) in enumerate(mesh.tess.edge_faces)
                              if len(e) == 2 and selen[i] > mnden
                              if sarea[e[0]] > mnden and sarea[e[1]] > mnden])    
      m       = len(e)
      (fis,q) = np.unique(np.concatenate([s,t]), return_inverse=True)
      (s,t)   = np.reshape(q, (2,-1))
      faces   = faces[fis]
      sarea   = sarea[fis]
      selen   = selen[e]
      (u,v)   = mesh.tess.indexed_edges[:,e]
      # we will use visual mag instead of cortical mag: this way we aren't worried about
      # varea going to 0 and creating a singularity, and it should have a linear 
      # relationship with eccentricity
      velen2  = (x[u] - x[v])**2 + (y[u] - y[v])**2
      vme     = velen2 / selen**2 # visual magnification: edges
      varea   = op.signed_face_areas(faces)
      vmf     = varea / sarea # visual magnification: faces
      vms     = vmf[s]
      vmt     = vmf[t]
      vsgns   = op.sign(vmf)
      f_magn  = (1.0 / m) * op.sum((vms - vsgns[t]*vme)**2 + (vmt - vsgns[s]*vme)**2)
      # [3] we want a special function for faces whose vmags are different signs
      if global_field_sign is None:
          f_sign = op.compose(op.piecewise(0, ((-np.inf, 0), -op.identity)), vms*vmt)
          f_sign = (1.0 / m) * op.sum(f_sign)
      else:
          vmfsgn = vmf * global_field_sign
          f_sign = op.compose(op.piecewise(0, ((-np.inf, 0), -op.identity)), vmfsgn)
          f_sign = (1.0 / m) * op.sum(f_sign)
      # and the edge potential...
      ex      = 0.5*(x[u] + x[v])
      ey      = 0.5*(y[u] + y[v])
      eecc2   = (ex**2 + ey**2)
      f_edge  = (1.0 / m) * op.sum(((x[u] - x[v])**2 + (y[u] - y[v])**2) / (eecc2 + 0.05))
      # This is the potential function:
      (k_meas, k_magn, k_sign, k_edge) = [
          0 if knob is None else (2**knob)
          for knob in (measurement_knob, magnification_knob, fieldsign_knob, edge_knob)]
      fs = (k_meas*f_meas, k_magn*f_magn, k_sign*f_sign, k_edge*f_edge)
      f = (fs[0] + fs[1] + fs[2] + fs[3]) / (k_meas + k_magn + k_sign + k_edge)
      xy0 = np.reshape(xy0, (-1,2))
      object.__setattr__(f, 'meta_data',
                         pyr.m(f_meas=f_meas, f_magn=f_magn, f_sign=f_sign, f_edge=f_edge,
                               mesh=mesh, X0=xy0))
      return f