def clean_retinotopy_potential(hemi, retinotopy=Ellipsis, mask=Ellipsis, weight=Ellipsis,<EOL>surface='<STR_LIT>', min_weight=Ellipsis, min_eccentricity=<NUM_LIT>,<EOL>visual_area=None, map_visual_areas=Ellipsis,<EOL>visual_area_field_signs=Ellipsis,<EOL>measurement_uncertainty=<NUM_LIT>, measurement_knob=<NUM_LIT:1>,<EOL>magnification_knob=<NUM_LIT:0>, fieldsign_knob=<NUM_LIT:6>, edge_knob=<NUM_LIT:0>):
from neuropythy.util import curry<EOL>import neuropythy.optimize as op<EOL>mesh = geo.to_mesh((hemi, surface))<EOL>rdat = (retinotopy_data(mesh) if retinotopy is Ellipsis   else<EOL>retinotopy            if pimms.is_map(retinotopy) else<EOL>retinotopy_data(mesh, retinotopy))<EOL>lbls = (rdat.get('<STR_LIT>') if visual_area is Ellipsis else<EOL>None                    if visual_area is None     else<EOL>hemi.property(visual_area))<EOL>wght = (rdat.get('<STR_LIT>') if weight is Ellipsis                      else<EOL>weight                         if pimms.is_vector(weight)                 else<EOL>rdat.get(weight)               if pimms.is_str(weight) and weight in rdat else<EOL>hemi.property(weight)          if weight is not None                      else<EOL>None)<EOL>if mask is Ellipsis:<EOL><INDENT>if mesh.coordinates.shape[<NUM_LIT:0>] == <NUM_LIT:2>: mask = mesh.indices<EOL>else: mask = geo.to_flatmap('<STR_LIT>', hemi, radius=np.pi/<NUM_LIT>).labels<EOL><DEDENT>else: mask = hemi.mask(mask, indices=True)<EOL>global_field_sign = None<EOL>if map_visual_areas and lbls is not None:<EOL><INDENT>vas = (np.unique(lbls)                    if map_visual_areas == '<STR_LIT:all>'           else<EOL>np.setdiff1d(np.unique(lbls), [<NUM_LIT:0>]) if map_visual_areas in [True,Ellipsis] else<EOL>np.unique(map_visual_areas))<EOL>if visual_area_field_signs is None: visual_area_field_signs = {}<EOL>elif visual_area_field_signs is Ellipsis: visual_area_field_signs = {<NUM_LIT:1>:-<NUM_LIT:1>, <NUM_LIT:2>:<NUM_LIT:1>, <NUM_LIT:3>:-<NUM_LIT:1>, <NUM_LIT:4>:<NUM_LIT:1>}<EOL>if pimms.is_int(map_visual_areas) or pimms.is_str(map_visual_areas):<EOL><INDENT>mask = np.intersect1d(mask, np.where(lbls == map_visual_areas)[<NUM_LIT:0>])<EOL>global_field_sign = visual_area_field_signs.get(map_visual_areas)<EOL><DEDENT>else: <EOL><INDENT>kw = dict(retinotopy=rdat, mask=mask, weight=wght,<EOL>surface=surface, min_weight=min_weight, min_eccentricity=min_eccentricity,<EOL>visual_area=lbls, measurement_uncertainty=measurement_uncertainty,<EOL>measurement_knob=measurement_knob,<EOL>magnification_knob=magnification_knob, fieldsign_knob=fieldsign_knob,<EOL>edge_knob=edge_knob, visual_area_field_signs=visual_area_field_signs)<EOL>return pimms.lazy_map(<EOL>{lbl: curry(clean_retinotopy_potential, hemi, map_visual_areas=lbl, **kw)<EOL>for lbl in vas})<EOL><DEDENT><DEDENT>(supermesh, orig_mask) = (mesh, mask)<EOL>rdat = {k:(v[mask] if len(v) > len(mask) else v) for (k,v) in six.iteritems(rdat)}<EOL>mesh = supermesh.submesh(mask)<EOL>mask = supermesh.tess.index(mesh.labels)<EOL>if len(mask) == len(orig_mask): smsk = np.arange(len(mask))<EOL>else:<EOL><INDENT>tmp = set(mask)<EOL>smsk = np.asarray([k for (k,u) in enumerate(orig_mask) if u in tmp])<EOL><DEDENT>n = mesh.vertex_count <EOL>N = <NUM_LIT:2>*n <EOL>if wght is None:                  wght = np.ones(n)<EOL>elif len(wght) == len(orig_mask): wght = np.array(wght)[smsk]<EOL>elif len(wght) > n:               wght = np.array(wght)[mask]<EOL>else:                             wght = np.array(wght)<EOL>wght[~np.isfinite(wght)] = <NUM_LIT:0><EOL>if min_eccentricity < <NUM_LIT:0> or np.isclose(min_eccentricity, <NUM_LIT:0>):<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>xy      = op.identity<EOL>(x,y)   = (xy[np.arange(<NUM_LIT:0>,N,<NUM_LIT:2>)],xy[np.arange(<NUM_LIT:1>,N,<NUM_LIT:2>)])<EOL>xy0     = np.array(as_retinotopy(rdat, '<STR_LIT>')).T<EOL>if len(xy0) == len(orig_mask): xy0 = xy0[smsk]<EOL>elif len(xy0) > n:             xy0 = xy0[mask]<EOL>(x0,y0) = xy0.T<EOL>xy0     = xy0.flatten()<EOL>ecc0    = np.sqrt(x0**<NUM_LIT:2> + y0**<NUM_LIT:2>)<EOL>ii      = np.where(ecc0 > min_eccentricity)[<NUM_LIT:0>]<EOL>minw    = (<NUM_LIT:0>                             if min_weight is None          else<EOL>np.percentile(wght[ii], <NUM_LIT:5>)    if min_weight is Ellipsis      else<EOL>min_weight                    if pimms.is_number(min_weight) else<EOL><NUM_LIT:0>                             if np.std(wght[ii]) < <NUM_LIT>  else<EOL>np.percentile(wght[ii], min_weight[<NUM_LIT:0>]))<EOL>ii      = np.intersect1d(ii, np.where(wght > minw)[<NUM_LIT:0>])<EOL>wsum    = np.sum(wght[ii])<EOL>if wsum < <NUM_LIT:0> or np.isclose(wsum, <NUM_LIT:0>): raise ValueError('<STR_LIT>')<EOL>wght    = wght / wsum<EOL>s2_meas = (measurement_uncertainty * ecc0[ii])**<NUM_LIT:2><EOL>d2_meas = (x[ii] - x0[ii])**<NUM_LIT:2> + (y[ii] - y0[ii])**<NUM_LIT:2><EOL>f_meas  = (<NUM_LIT:1> - op.exp(-<NUM_LIT:0.5>*d2_meas/s2_meas)) * wght[ii]<EOL>f_meas  = op.sum(f_meas)<EOL>sarea   = mesh.face_areas<EOL>faces   = mesh.tess.indexed_faces.T<EOL>selen   = mesh.edge_lengths<EOL>mnden   = <NUM_LIT><EOL>(e,s,t) = np.transpose([(i,e[<NUM_LIT:0>],e[<NUM_LIT:1>]) for (i,e) in enumerate(mesh.tess.edge_faces)<EOL>if len(e) == <NUM_LIT:2> and selen[i] > mnden<EOL>if sarea[e[<NUM_LIT:0>]] > mnden and sarea[e[<NUM_LIT:1>]] > mnden])    <EOL>m       = len(e)<EOL>(fis,q) = np.unique(np.concatenate([s,t]), return_inverse=True)<EOL>(s,t)   = np.reshape(q, (<NUM_LIT:2>,-<NUM_LIT:1>))<EOL>faces   = faces[fis]<EOL>sarea   = sarea[fis]<EOL>selen   = selen[e]<EOL>(u,v)   = mesh.tess.indexed_edges[:,e]<EOL>velen2  = (x[u] - x[v])**<NUM_LIT:2> + (y[u] - y[v])**<NUM_LIT:2><EOL>vme     = velen2 / selen**<NUM_LIT:2> <EOL>varea   = op.signed_face_areas(faces)<EOL>vmf     = varea / sarea <EOL>vms     = vmf[s]<EOL>vmt     = vmf[t]<EOL>vsgns   = op.sign(vmf)<EOL>f_magn  = (<NUM_LIT:1.0> / m) * op.sum((vms - vsgns[t]*vme)**<NUM_LIT:2> + (vmt - vsgns[s]*vme)**<NUM_LIT:2>)<EOL>if global_field_sign is None:<EOL><INDENT>f_sign = op.compose(op.piecewise(<NUM_LIT:0>, ((-np.inf, <NUM_LIT:0>), -op.identity)), vms*vmt)<EOL>f_sign = (<NUM_LIT:1.0> / m) * op.sum(f_sign)<EOL><DEDENT>else:<EOL><INDENT>vmfsgn = vmf * global_field_sign<EOL>f_sign = op.compose(op.piecewise(<NUM_LIT:0>, ((-np.inf, <NUM_LIT:0>), -op.identity)), vmfsgn)<EOL>f_sign = (<NUM_LIT:1.0> / m) * op.sum(f_sign)<EOL><DEDENT>ex      = <NUM_LIT:0.5>*(x[u] + x[v])<EOL>ey      = <NUM_LIT:0.5>*(y[u] + y[v])<EOL>eecc2   = (ex**<NUM_LIT:2> + ey**<NUM_LIT:2>)<EOL>f_edge  = (<NUM_LIT:1.0> / m) * op.sum(((x[u] - x[v])**<NUM_LIT:2> + (y[u] - y[v])**<NUM_LIT:2>) / (eecc2 + <NUM_LIT>))<EOL>(k_meas, k_magn, k_sign, k_edge) = [<EOL><NUM_LIT:0> if knob is None else (<NUM_LIT:2>**knob)<EOL>for knob in (measurement_knob, magnification_knob, fieldsign_knob, edge_knob)]<EOL>fs = (k_meas*f_meas, k_magn*f_magn, k_sign*f_sign, k_edge*f_edge)<EOL>f = (fs[<NUM_LIT:0>] + fs[<NUM_LIT:1>] + fs[<NUM_LIT:2>] + fs[<NUM_LIT:3>]) / (k_meas + k_magn + k_sign + k_edge)<EOL>xy0 = np.reshape(xy0, (-<NUM_LIT:1>,<NUM_LIT:2>))<EOL>object.__setattr__(f, '<STR_LIT>',<EOL>pyr.m(f_meas=f_meas, f_magn=f_magn, f_sign=f_sign, f_edge=f_edge,<EOL>mesh=mesh, X0=xy0))<EOL>return f<EOL>
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.