def f2():
    D = np.zeros((3, 3))
    D[:] = np.nan

    A, B, C = 0, 1, 2

    ab, bc, ca = 3, 4, 5

    names = {A: 'A', B: 'B', C: 'C', ab: 'ab', bc: 'bc', ca: 'ca'}

    d = np.pi / 25
    D[A, B] = D[B, C] = D[C, A] = d

    fillD(D)

#    some = svds(np.cos(D), 10)
#    
    S = best_embedding_on_sphere(np.cos(D), 3)

    S2 = np.zeros((3, 6))
    for i in range(3):
        S2[:, i] = S[:, i]

    S2[:, ab] = midpoint_on_sphere(S[:, A], S[:, B])
    S2[:, bc] = midpoint_on_sphere(S[:, B], S[:, C])
    S2[:, ca] = midpoint_on_sphere(S[:, C], S[:, A])

    D2 = distances_from_directions(S2)

    @contract(x='array[MxN]')
    def pprint(x, format='%.3f'): #@ReservedAssignment
        for a in range(x.shape[0]):
            sys.stdout.write('[')
            for b in range(x.shape[1]):
                sys.stdout.write(format % x[a, b])
                sys.stdout.write(', ')
            sys.stdout.write(']\n')
    pprint(D2, '%.2f')

    S3 = best_embedding_on_sphere(np.cos(D2), 3)
    D2 = distances_from_directions(S3)

    assert_allclose(D[A, B], 2 * D2[A, ab])
    assert_allclose(D[A, C], 2 * D2[A, ca])
    assert_allclose(D[C, B], 2 * D2[C, bc])

    for t in [[A, B, C],
              [ab, bc, ca],
              [A, ab, ca],
              [B, ab, bc],
              [C, ca, bc]]:
        x, y, z = t
        d = [D2[x, y], D2[y, z], D2[z, x]]
        equil = allclose(d[0], d[1]) and allclose(d[1], d[2])
        name = '-'.join(names[x] for x in t)
        mark = '(equiv)' if equil else ""
        print('Triangle %s: %.2f %.2f %.2f %s'
              % (name, d[0], d[1], d[2], mark))

 def is_state_dist(self, state_dist):
     for s in state_dist:
         self.is_state(s)
     p = sum(state_dist.values())
     p = float(p)
     if not allclose(p, 1.0):
         raise ValueError('PD sums to %f' % p)

 def modify_avg_degree(self, value):
     """
     Modifies (increases) average degree based on given value by
     modifying nodes commRange."""
     # assert all nodes have same commRange
     assert allclose([n.commRange for n in self], self.nodes()[0].commRange)
     #TODO: implement decreasing of degree, preserve connected network
     assert value + settings.DEG_ATOL > self.avg_degree()  # only increment
     step_factor = 7.
     steps = [0]
     #TODO: while condition should call validate
     while not allclose(self.avg_degree(), value, atol=settings.DEG_ATOL):
         steps.append((value - self.avg_degree())*step_factor)
         for node in self:
             node.commRange += steps[-1]
         # variable step_factor for step size for over/undershoot cases
         if len(steps)>2 and sign(steps[-2])!=sign(steps[-1]):
             step_factor /= 2
     logger.debug("Modified degree to %f" % self.avg_degree())

def polydiv(u, v):
    """
    Returns the quotient and remainder of polynomial division.

    The input arrays specify the polynomial terms in turn with a length equal
    to the polynomial degree plus 1.

    Parameters
    ----------
    u : {array_like, poly1d}
        Dividend polynomial.
    v : {array_like, poly1d}
        Divisor polynomial.

    Returns
    -------
    q : ndarray
        Polynomial terms of quotient.
    r : ndarray
        Remainder of polynomial division.

    See Also
    --------
    poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub,
    polyval

    Examples
    --------
    .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25

    >>> x = np.array([3.0, 5.0, 2.0])
    >>> y = np.array([2.0, 1.0])
    >>> np.polydiv(x, y)
    >>> (array([ 1.5 ,  1.75]), array([ 0.25]))

    """
    truepoly = (isinstance(u, poly1d) or isinstance(u, poly1d))
    u = atleast_1d(u) + 0.0
    v = atleast_1d(v) + 0.0
    # w has the common type
    w = u[0] + v[0]
    m = len(u) - 1
    n = len(v) - 1
    scale = 1. / v[0]
    q = NX.zeros((max(m - n + 1, 1),), w.dtype)
    r = u.copy()
    for k in range(0, m-n+1):
        d = scale * r[k]
        q[k] = d
        r[k:k+n+1] -= d*v
    while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
        r = r[1:]
    if truepoly:
        return poly1d(q), poly1d(r)
    return q, r

def real_if_close(a, tol=100):
    """
    If complex input returns a real array if complex parts are close to zero.

    "Close to zero" is defined as `tol` * (machine epsilon of the type for
    `a`).

    Parameters
    ----------
    a : array_like
        Input array.
    tol : float
        Tolerance in machine epsilons for the complex part of the elements
        in the array.

    Returns
    -------
    out : ndarray
        If `a` is real, the type of `a` is used for the output.  If `a`
        has complex elements, the returned type is float.

    See Also
    --------
    real, imag, angle

    Notes
    -----
    Machine epsilon varies from machine to machine and between data types
    but Python floats on most platforms have a machine epsilon equal to
    2.2204460492503131e-16.  You can use 'np.finfo(np.float).eps' to print
    out the machine epsilon for floats.

    Examples
    --------
    >>> np.finfo(np.float).eps
    2.2204460492503131e-16

    >>> np.real_if_close([2.1 + 4e-14j], tol=1000)
    array([ 2.1])
    >>> np.real_if_close([2.1 + 4e-13j], tol=1000)
    array([ 2.1 +4.00000000e-13j])

    """
    a = asanyarray(a)
    if not issubclass(a.dtype.type, _nx.complexfloating):
        return a
    if tol > 1:
        from numpy.core import getlimits

        f = getlimits.finfo(a.dtype.type)
        tol = f.eps * tol
    if _nx.allclose(a.imag, 0, atol=tol):
        a = a.real
    return a

def plot_vertical_colorbar(pl, colorbar, bar_x, bar_w, bar_y_min, bar_y_max,
                           vdist, min_value, max_value, scale_format=None):
    
    if scale_format is None:
        if isinstance(max_value, int):
            scale_format = '%d'
        else:
            scale_format = '%.2f'
            
        if allclose(max_value, +1) and allclose(min_value, -1):
            scale_format = '%+d'
        
    label_min = scale_format % min_value
    if min_value == 0:
        label_min = '0'
    label_max = scale_format % max_value
                                   
    ex = [bar_x - bar_w, bar_x + bar_w, bar_y_min, bar_y_max]

    border = colorbar.shape[1] / 10
    print('using border %d for %s' % (border, colorbar.shape))
    for b in range(border):
        colorbar[0 + b, :, 0:3] = 0
        colorbar[-1 - b, :, 0:3] = 0
        colorbar[:, 0 + b, 0:3] = 0
        colorbar[:, -1 - b, 0:3] = 0
         
    x = pl.imshow(colorbar, origin='lower',
              extent=ex, aspect='auto',
              interpolation='nearest')
    
    x.set_clip_on(False)
#    pl.fill([ex[0], ex[0], ex[1], ex[1]],
#            [ex[2], ex[3], ex[3], ex[2]], facecolor='none', edgecolor='k')

    pl.annotate(label_min, (bar_x, bar_y_min - vdist),
                horizontalalignment='center', verticalalignment='top')
    pl.annotate(label_max, (bar_x, bar_y_max + vdist),
                horizontalalignment='center', verticalalignment='bottom') 

def real_if_close(a,tol=100):
    """If a is a complex array, return it as a real array if the imaginary
    part is close enough to zero.

    "Close enough" is defined as tol*(machine epsilon of a's element type).
    """
    a = asanyarray(a)
    if not issubclass(a.dtype.type, _nx.complexfloating):
        return a
    if tol > 1:
        import getlimits
        f = getlimits.finfo(a.dtype.type)
        tol = f.eps * tol
    if _nx.allclose(a.imag, 0, atol=tol):
        a = a.real
    return a

 def validate_params(self, params):
     """ Validate if given network params match its real params. """
     logger.info('Validating params')
     count = params.get('count', None)  #  for unit tests
     if count:
         if isinstance(count, list):
             assert(len(self) in count)
         else:
             assert(len(self)==count)
     n_min = params.get('n_min', 0)
     n_max = params.get('n_max', Inf)
     assert(len(self)>=n_min and len(self)<=n_max)
     for param, value in params.items():
         if param=='connected':
             assert(not value or is_connected(self))
         elif param=='degree':
             assert(allclose(self.avg_degree(), value,
                             atol=settings.DEG_ATOL))
         elif param=='environment':
             assert(self.environment.__class__==value.__class__)
         elif param=='channelType':
             assert(self.channelType.__class__==value.__class__)
         elif param=='comm_range':
             for node in self:
                 assert(node.commRange==value)
         elif param=='sensors':
             compositeSensor = CompositeSensor(Node(), value)
             for node in self:
                 assert(all(map(lambda s1, s2: pymote_equal_objects(s1, s2),
                                node.sensors, compositeSensor.sensors)))
         elif param=='aoa_pf_scale':
             for node in self:
                 for sensor in node.sensors:
                     if sensor.name()=='AoASensor':
                         assert(sensor.probabilityFunction.scale==value)
         elif param=='dist_pf_scale':
             for node in self:
                 for sensor in node.sensors:
                     if sensor.name()=='DistSensor':
                         assert(sensor.probabilityFunction.scale==value)
         #TODO: refactor this part as setting algorithms resets nodes
         """

 def validate_params(self, params):
     """ Validate if given network params match its real params. """
     logger.info("Validating params")
     count = params.get("count", None)  #  for unit tests
     if count:
         if isinstance(count, list):
             assert len(self) in count
         else:
             assert len(self) == count
     n_min = params.get("n_min", 0)
     n_max = params.get("n_max", Inf)
     assert len(self) >= n_min and len(self) <= n_max
     for param, value in params.items():
         if param == "connected":
             assert not value or is_connected(self)
         elif param == "degree":
             assert allclose(self.avg_degree(), value, atol=settings.DEG_ATOL)
         elif param == "environment":
             assert self.environment.__class__ == value.__class__
         elif param == "channelType":
             assert self.channelType.__class__ == value.__class__
         elif param == "comm_range":
             for node in self:
                 assert node.commRange == value
         elif param == "sensors":
             compositeSensor = CompositeSensor(Node(), value)
             for node in self:
                 assert all(map(lambda s1, s2: pymote_equal_objects(s1, s2), node.sensors, compositeSensor.sensors))
         elif param == "aoa_pf_scale":
             for node in self:
                 for sensor in node.sensors:
                     if sensor.name() == "AoASensor":
                         assert sensor.probabilityFunction.scale == value
         elif param == "dist_pf_scale":
             for node in self:
                 for sensor in node.sensors:
                     if sensor.name() == "DistSensor":
                         assert sensor.probabilityFunction.scale == value
         # TODO: refactor this part as setting algorithms resets nodes
         """

def polydiv(u, v):
    """Computes q and r polynomials so that u(s) = q(s)*v(s) + r(s)
    and deg r < deg v.
    """
    truepoly = (isinstance(u, poly1d) or isinstance(u, poly1d))
    u = atleast_1d(u)
    v = atleast_1d(v)
    m = len(u) - 1
    n = len(v) - 1
    scale = 1. / v[0]
    q = NX.zeros((max(m-n+1,1),), float)
    r = u.copy()
    for k in range(0, m-n+1):
        d = scale * r[k]
        q[k] = d
        r[k:k+n+1] -= d*v
    while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
        r = r[1:]
    if truepoly:
        q = poly1d(q)
        r = poly1d(r)
    return q, r

def find_polytope_bounds_after_linear(streamels, A, streamels2):
    ''' 
        Fills the 'lower' and 'upper' part of streamels2
        given the A transform.
    '''
    M, N = A.shape
    check_streamels_1D_size(streamels, N)
    check_streamels_1D_size(streamels2, M)

    is_diagonal = (M == N) and allclose(np.diag(np.diagonal(A)), A)

    if is_diagonal:
        streamels2['lower'] = np.dot(A, streamels['lower'])
        streamels2['upper'] = np.dot(A, streamels['upper'])

    else:  # TODO: do something smarter here
        norm = np.abs(A).sum()  # XXX
        lb = np.max(np.abs(streamels['lower']))
        ub = np.max(np.abs(streamels['upper']))
        B = max(lb, ub)
        streamels2['lower'] = -B * norm
        streamels2['upper'] = +B * norm

def polydiv(u, v):
    """
    Returns the quotient and remainder of polynomial division.

    The input arrays are the coefficients (including any coefficients
    equal to zero) of the "numerator" (dividend) and "denominator"
    (divisor) polynomials, respectively.

    Parameters
    ----------
    u : array_like or poly1d
        Dividend polynomial's coefficients.

    v : array_like or poly1d
        Divisor polynomial's coefficients.

    Returns
    -------
    q : ndarray
        Coefficients, including those equal to zero, of the quotient.
    r : ndarray
        Coefficients, including those equal to zero, of the remainder.

    See Also
    --------
    poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub,
    polyval

    Notes
    -----
    Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need
    not equal `v.ndim`. In other words, all four possible combinations -
    ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``,
    ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work.

    Examples
    --------
    .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25

    >>> x = np.array([3.0, 5.0, 2.0])
    >>> y = np.array([2.0, 1.0])
    >>> np.polydiv(x, y)
    (array([ 1.5 ,  1.75]), array([ 0.25]))

    """
    truepoly = (isinstance(u, poly1d) or isinstance(u, poly1d))
    u = atleast_1d(u) + 0.0
    v = atleast_1d(v) + 0.0
    # w has the common type
    w = u[0] + v[0]
    m = len(u) - 1
    n = len(v) - 1
    scale = 1. / v[0]
    q = NX.zeros((max(m - n + 1, 1),), w.dtype)
    r = u.copy()
    for k in range(0, m-n+1):
        d = scale * r[k]
        q[k] = d
        r[k:k+n+1] -= d*v
    while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
        r = r[1:]
    if truepoly:
        return poly1d(q), poly1d(r)
    return q, r

def midpoint_on_sphere(s1, s2):
    if allclose(s1, -s2): # TODO: add precision
        raise ValueError()
    else:
        v = (s1 + s2) * 0.5
        return v / np.linalg.norm(v)