def retinotopy_anchors(mesh, mdl,
                       polar_angle=None, eccentricity=None,
                       weight=None, weight_min=0.1,
                       field_sign_weight=0, field_sign=None, invert_field_sign=False,
                       radius_weight=0, radius_weight_source='Wandell2015', radius=None,
                       model_field_sign=None,
                       model_hemi=Ellipsis,
                       scale=1,
                       shape='Gaussian', suffix=None,
                       sigma=[0.1, 2.0, 8.0],
                       select='close'):
    '''
    retinotopy_anchors(mesh, model) is intended for use with the mesh_register function and the
    retinotopy_model() function and/or the RetinotopyModel class; it yields a description of the
    anchor points that tie relevant vertices the given mesh to points predicted by the given model
    object. Any instance of the RetinotopyModel class should work as a model argument; this includes
    SchiraModel objects as well as RetinotopyMeshModel objects such as those returned by the
    retinotopy_model() function. If the model given is a string, then it is passed to the
    retinotopy_model() function first.

    Options:
      * polar_angle (default None) specifies that the given data should be used in place of the
        'polar_angle' or 'PRF_polar_angle'  property values. The given argument must be numeric and
        the same length as the the number of vertices in the mesh. If None is given, then the
        property value of the mesh is used; if a list is given and any element is None, then the
        weight for that vertex is treated as a zero. If the option is a string, then the property
        value with the same name isused as the polar_angle data.
      * eccentricity (default None) specifies that the given data should be used in places of the
        'eccentricity' or 'PRF_eccentricity' property values. The eccentricity option is handled 
        virtually identically to the polar_angle option.
      * weight (default None) specifies that the weight or scale of the data; this is handled
        generally like the polar_angle and eccentricity options, but may also be 1, indicating that
        all vertices with polar_angle and eccentricity values defined will be given a weight of 1.
        If weight is left as None, then the function will check for 'weight',
        'variance_explained', 'PRF_variance_explained', and 'retinotopy_weight' values and will use
        the first found (in that order). If none of these is found, then a value of 1 is assumed.
      * weight_min (default 0) specifies that the weight must be higher than the given value inn
        order to be included in the fit; vertices with weights below this value have their weights
        truncated to 0.
      * scale (default 1) specifies a constant by which to multiply all weights for all anchors; the
        value None is interpreted as 1.
      * shape (default 'Gaussian') specifies the shape of the potential function (see mesh_register)
      * model_hemi (default: None) specifies the hemisphere of the model to load; if None, then
        looks for a non-specific model.
      * suffix (default None) specifies any additional arguments that should be appended to the 
        potential function description list that is produced by this function; i.e., the 
        retinotopy_anchors function produces a list, and the contents of suffix, if given and not
        None, are appended to that list (see mesh_register).
      * select (default 'close') specifies a function that will be called with two arguments for
        every vertex given an anchor; the arguments are the vertex label and the matrix of anchors.
        The function should return a list of anchors to use for the label (None is equivalent to
        lambda id,anc: anc). The parameter may alternately be specified using the string 'close':
        select=['close', [k]] indicates that any anchor more than k times the average edge-length in
        the mesh should be excluded; a value of just ['close', k] on the other hand indicates that
        any anchor more than k distance from the vertex should be exlcuded. The default value,
        'close', is equivalent to ['close', [40]].
      * sigma (default [0.1, 2.0, 4.0]) specifies how the sigma parameter should be handled; if
        None, then no sigma value is specified; if a single number, then all sigma values are
        assigned that value; if a list of three numbers, then the first is the minimum sigma value,
        the second is the fraction of the minimum distance between paired anchor points, and the 
        last is the maximum sigma --- the idea with this form of the argument is that the ideal
        sigma value in many cases is approximately 0.25 to 0.5 times the distance between anchors
        to which a single vertex is attracted; for any anchor a to which a vertex u is attracted,
        the sigma of a is the middle sigma-argument value times the minimum distance from a to all
        other anchors to which u is attracted (clipped by the min and max sigma).
      * field_sign_weight (default: 0) specifies the amount of weight that should be put on the
        retinotopic field of the model as a method of attenuating the weights on those anchors whose
        empirical retinotopic values and predicted model locations do not match. The weight that
        results is calculated from the difference in empirical field-sign for each vertex and the
        visual area field sign based on the labels in the model. The higher the field-sign weight,
        (approaching 1) the more the resulting value is a geometric mean of the field-sign-based
        weight and the original weights. As this value approaches 0, the resulting weights are more
        like the original weights.
      * radius_weight (default: 0) specifies the amount of weight that should be put on the
        receptive field radius of the model as a method of attenuating the weights on those anchors
        whose empirical retinotopic values and predicted model locations do not match. The weight
        that results is calculated from the difference in empirical RF radius for each vertex and
        the predicted RF radius based on the labels in the model. The higher the radius weight,
        (approaching 1) the more the resulting value is a geometric mean of the field-sign-based
        weight and the original weights. As this value approaches 0, the resulting weights are more
        like the original weights.
      * radius_weight_source (default: 'Wandell2015') specifies the source for predicting RF radius;
        based on eccentricity and visual area label.

    Example:
     # The retinotopy_anchors function is intended for use with mesh_register, as follows:
     # Define our Schira Model:
     model = neuropythy.registration.SchiraModel()
     # Make sure our mesh has polar angle, eccentricity, and weight data:
     mesh.prop('polar_angle',  polar_angle_vertex_data);
     mesh.prop('eccentricity', eccentricity_vertex_data);
     mesh.prop('weight',       variance_explained_vertex_data);
     # register the mesh using the retinotopy and model:
     registered_mesh = neuropythy.registration.mesh_register(
        mesh,
        ['mesh', retinotopy_anchors(mesh, model)],
        max_step_size=0.05,
        max_steps=2000)
    '''
    if pimms.is_str(mdl):
        hemi = None
        if pimms.is_str(model_hemi):
            model_hemi = model_hemi.upper()
            hemnames = {k:h
                        for (h,als) in [('LH', ['LH','L','LEFT','RHX','RX']),
                                        ('RH', ['RH','R','RIGHT','LHX','LX'])]
                        for k in als}
            if model_hemi in hemnames: hemi = hemnames[model_hemi]
            else: raise ValueError('Unrecognized hemisphere name: %s' % model_hemi)
        elif model_hemi is not None:
            raise ValueError('model_hemi must be a string, Ellipsis, or None')
        mdl = retinotopy_model(mdl, hemi=hemi)
    if not isinstance(mdl, RetinotopyModel):
        raise RuntimeError('given model is not a RetinotopyModel instance!')
    if not isinstance(mesh, geo.Mesh):
        raise RuntimeError('given mesh is not a Mesh object!')
    n = mesh.vertex_count
    X = mesh.coordinates.T
    if weight_min is None: weight_min = 0
    # make sure we have our polar angle/eccen/weight values:
    # (weight is odd because it might be a single number, so handle that first)
    (polar_angle, eccentricity, weight) = [
        extract_retinotopy_argument(mesh, name, arg, default='empirical')
        for (name, arg) in [
                ('polar_angle', polar_angle),
                ('eccentricity', eccentricity),
                ('weight', np.full(n, weight) if pimms.is_number(weight) else weight)]]
    # Make sure they contain no None/invalid values
    (polar_angle, eccentricity, weight) = _retinotopy_vectors_to_float(
        polar_angle, eccentricity, weight,
        weight_min=weight_min)
    if np.sum(weight > 0) == 0:
        raise ValueError('No positive weights found')
    idcs = np.where(weight > 0)[0]
    # Interpret the select arg if necessary (but don't apply it yet)
    select = ['close', [40]] if select == 'close'   else \
             ['close', [40]] if select == ['close'] else \
             select
    if select is None:
        select = lambda a,b: b
    elif ((pimms.is_vector(select) or is_list(select) or is_tuple(select))
          and len(select) == 2 and select[0] == 'close'):
        if pimms.is_vector(select[1]): d = np.mean(mesh.edge_lengths) * select[1][0]
        else:                          d = select[1]
        select = lambda idx,ancs: [a for a in ancs if a[0] is not None if npla.norm(X[idx] - a) < d]
    # Okay, apply the model:
    res = mdl.angle_to_cortex(polar_angle[idcs], eccentricity[idcs])
    oks = np.isfinite(np.sum(np.reshape(res, (res.shape[0], -1)), axis=1))
    # Organize the data; trim out those not selected
    data = [[[i for _ in r], r, [ksidx[tuple(a)] for a in r]]
            for (i,r0,ok) in zip(idcs, res, oks) if ok
            for ksidx in [{tuple(a):(k+1) for (k,a) in enumerate(r0)}]
            for r in [select(i, r0)]
            if len(r) > 0]
    # Flatten out the data into arguments for Java
    idcs = [int(i) for d in data for i in d[0]]
    ancs = np.asarray([pt for d in data for pt in d[1]]).T
    labs = np.asarray([ii for d in data for ii in d[2]]).T
    # Get just the relevant weights and the scale
    wgts = np.asarray(weight[idcs] * (1 if scale is None else scale))
    # add in the field-sign weights and radius weights if requested here;
    if not np.isclose(field_sign_weight, 0) and mdl.area_name_to_id is not None:
        id2n = mdl.area_id_to_name
        if field_sign is True or field_sign is Ellipsis or field_sign is None:
            from .cmag import cmag
            r = {'polar_angle': polar_angle,  'eccentricity': eccentricity}
            #field_sign = retinotopic_field_sign(mesh, retinotopy=r)
            field_sign = cmag(mesh, r)['field_sign']
        elif pimms.is_str(field_sign): field_sign = mesh.prop(field_sign)
        field_sign = np.asarray(field_sign)
        if invert_field_sign: field_sign = -field_sign
        fswgts = 1.0 - 0.25 * np.asarray(
            [(fs - visual_area_field_signs[id2n[l]]) if l in id2n else 0
             for (l,fs) in zip(labs,field_sign[idcs])])**2
        # average the weights at some fraction with the original weights
        fswgts = field_sign_weight*fswgts + (1 - field_sign_weight)*wgts
    else: fswgts = None
    # add in radius weights if requested as well
    if not np.isclose(radius_weight, 0) and mdl.area_name_to_id is not None:
        id2n = mdl.area_id_to_name
        emprad = extract_retinotopy_argument(mesh, 'radius', radius, default='empirical')
        emprad = emprad[idcs]
        emprad = np.argsort(np.argsort(emprad)) * (1.0 / len(emprad)) - 0.5
        eccs = eccentricity[idcs]
        prerad = np.asarray([predict_pRF_radius(ecc, id2n[lbl], source=radius_weight_source)
                             for (ecc,lbl) in zip(eccs,labs)])
        prerad = np.argsort(np.argsort(prerad)) * (1.0 / len(prerad)) - 0.5
        rdwgts = 1.0 - (emprad - prerad)**2
        # average the weights at some fraction with the original weights
        rdwgts = radius_weight*rdwgts + (1-radius_weight)*wgts
    else: rdwgts = None
    # apply the weights
    if fswgts is not None:
        if rdwgts is not None: wgts = np.power(fswgts*rdwgts*wgts, 1.0/3.0)
        else:                  wgts = np.sqrt(fswgts*wgts)
    elif rdwgts is not None:   wgts = np.sqrt(rdwgts*wgts)
    # Figure out the sigma parameter:
    if sigma is None: sigs = None
    elif pimms.is_number(sigma): sigs = sigma
    elif pimms.is_vector(sigma) and len(sigma) == 3:
        [minsig, mult, maxsig] = sigma
        sigs = np.clip(
            [mult*min([npla.norm(a0 - a) for a in anchs if a is not a0]) if len(iii) > 1 else maxsig
             for (iii,anchs,_) in data
             for a0 in anchs],
            minsig, maxsig)
    else:
        raise ValueError('sigma must be a number or a list of 3 numbers')
    # okay, we've partially parsed the data that was given; now we can construct the final list of
    # instructions:
    tmp =  (['anchor', shape,
             np.asarray(idcs, dtype=np.int),
             np.asarray(ancs, dtype=np.float64),
             'scale', np.asarray(wgts, dtype=np.float64)]
            + ([] if sigs is None else ['sigma', sigs])
            + ([] if suffix is None else suffix))
    return tmp
retinotopy_anchors(mesh, model) is intended for use with the mesh_register function and the
    retinotopy_model() function and/or the RetinotopyModel class; it yields a description of the
    anchor points that tie relevant vertices the given mesh to points predicted by the given model
    object. Any instance of the RetinotopyModel class should work as a model argument; this includes
    SchiraModel objects as well as RetinotopyMeshModel objects such as those returned by the
    retinotopy_model() function. If the model given is a string, then it is passed to the
    retinotopy_model() function first.

    Options:
      * polar_angle (default None) specifies that the given data should be used in place of the
        'polar_angle' or 'PRF_polar_angle'  property values. The given argument must be numeric and
        the same length as the the number of vertices in the mesh. If None is given, then the
        property value of the mesh is used; if a list is given and any element is None, then the
        weight for that vertex is treated as a zero. If the option is a string, then the property
        value with the same name isused as the polar_angle data.
      * eccentricity (default None) specifies that the given data should be used in places of the
        'eccentricity' or 'PRF_eccentricity' property values. The eccentricity option is handled 
        virtually identically to the polar_angle option.
      * weight (default None) specifies that the weight or scale of the data; this is handled
        generally like the polar_angle and eccentricity options, but may also be 1, indicating that
        all vertices with polar_angle and eccentricity values defined will be given a weight of 1.
        If weight is left as None, then the function will check for 'weight',
        'variance_explained', 'PRF_variance_explained', and 'retinotopy_weight' values and will use
        the first found (in that order). If none of these is found, then a value of 1 is assumed.
      * weight_min (default 0) specifies that the weight must be higher than the given value inn
        order to be included in the fit; vertices with weights below this value have their weights
        truncated to 0.
      * scale (default 1) specifies a constant by which to multiply all weights for all anchors; the
        value None is interpreted as 1.
      * shape (default 'Gaussian') specifies the shape of the potential function (see mesh_register)
      * model_hemi (default: None) specifies the hemisphere of the model to load; if None, then
        looks for a non-specific model.
      * suffix (default None) specifies any additional arguments that should be appended to the 
        potential function description list that is produced by this function; i.e., the 
        retinotopy_anchors function produces a list, and the contents of suffix, if given and not
        None, are appended to that list (see mesh_register).
      * select (default 'close') specifies a function that will be called with two arguments for
        every vertex given an anchor; the arguments are the vertex label and the matrix of anchors.
        The function should return a list of anchors to use for the label (None is equivalent to
        lambda id,anc: anc). The parameter may alternately be specified using the string 'close':
        select=['close', [k]] indicates that any anchor more than k times the average edge-length in
        the mesh should be excluded; a value of just ['close', k] on the other hand indicates that
        any anchor more than k distance from the vertex should be exlcuded. The default value,
        'close', is equivalent to ['close', [40]].
      * sigma (default [0.1, 2.0, 4.0]) specifies how the sigma parameter should be handled; if
        None, then no sigma value is specified; if a single number, then all sigma values are
        assigned that value; if a list of three numbers, then the first is the minimum sigma value,
        the second is the fraction of the minimum distance between paired anchor points, and the 
        last is the maximum sigma --- the idea with this form of the argument is that the ideal
        sigma value in many cases is approximately 0.25 to 0.5 times the distance between anchors
        to which a single vertex is attracted; for any anchor a to which a vertex u is attracted,
        the sigma of a is the middle sigma-argument value times the minimum distance from a to all
        other anchors to which u is attracted (clipped by the min and max sigma).
      * field_sign_weight (default: 0) specifies the amount of weight that should be put on the
        retinotopic field of the model as a method of attenuating the weights on those anchors whose
        empirical retinotopic values and predicted model locations do not match. The weight that
        results is calculated from the difference in empirical field-sign for each vertex and the
        visual area field sign based on the labels in the model. The higher the field-sign weight,
        (approaching 1) the more the resulting value is a geometric mean of the field-sign-based
        weight and the original weights. As this value approaches 0, the resulting weights are more
        like the original weights.
      * radius_weight (default: 0) specifies the amount of weight that should be put on the
        receptive field radius of the model as a method of attenuating the weights on those anchors
        whose empirical retinotopic values and predicted model locations do not match. The weight
        that results is calculated from the difference in empirical RF radius for each vertex and
        the predicted RF radius based on the labels in the model. The higher the radius weight,
        (approaching 1) the more the resulting value is a geometric mean of the field-sign-based
        weight and the original weights. As this value approaches 0, the resulting weights are more
        like the original weights.
      * radius_weight_source (default: 'Wandell2015') specifies the source for predicting RF radius;
        based on eccentricity and visual area label.

    Example:
     # The retinotopy_anchors function is intended for use with mesh_register, as follows:
     # Define our Schira Model:
     model = neuropythy.registration.SchiraModel()
     # Make sure our mesh has polar angle, eccentricity, and weight data:
     mesh.prop('polar_angle',  polar_angle_vertex_data);
     mesh.prop('eccentricity', eccentricity_vertex_data);
     mesh.prop('weight',       variance_explained_vertex_data);
     # register the mesh using the retinotopy and model:
     registered_mesh = neuropythy.registration.mesh_register(
        mesh,
        ['mesh', retinotopy_anchors(mesh, model)],
        max_step_size=0.05,
        max_steps=2000)