Python numpy.linalg 模块，norm() 实例源码

def test_gmmmap_swap():
    static_dim = 2
    T = 10
    windows = _get_windows_set()[-1]

    np.random.seed(1234)
    src_mc = np.random.rand(T, static_dim * len(windows))
    tgt_mc = np.random.rand(T, static_dim * len(windows))

    # pseudo parallel data
    XY = np.concatenate((src_mc, tgt_mc), axis=-1)
    gmm = GaussianMixture(n_components=4)
    gmm.fit(XY)

    paramgen = MLPG(gmm, windows=windows, swap=False)
    swap_paramgen = MLPG(gmm, windows=windows, swap=True)

    mc_converted1 = paramgen.transform(src_mc)
    mc_converted2 = swap_paramgen.transform(tgt_mc)

    src_mc = src_mc[:, :static_dim]
    tgt_mc = tgt_mc[:, :static_dim]

    assert norm(tgt_mc - mc_converted1) < norm(src_mc - mc_converted1)
    assert norm(tgt_mc - mc_converted2) > norm(src_mc - mc_converted2)

def kNN_entity(self, vec, topk=10, method=0, self_vec_id=None):
        q = []
        for i in range(len(self.vec_e)):
            #skip self
            if self_vec_id != None and i == self_vec_id:
                continue
            if method == 1:
                dist = SP.distance.cosine(vec, self.vec_e[i])
            else:
                dist = LA.norm(vec - self.vec_e[i])
            if len(q) < topk:
                HP.heappush(q, self.index_dist(i, dist))
            else:
                #indeed it fetches the biggest
                tmp = HP.nsmallest(1, q)[0]
                if tmp.dist > dist:
                    HP.heapreplace(q, self.index_dist(i, dist) )
        rst = []
        while len(q) > 0:
            item = HP.heappop(q)
            rst.insert(0, (self.vocab_e[self.vec2e[item.index]], item.dist))
        return rst

    #given entity name, find kNN

def _ncc_c(x, y):
    """
    >>> _ncc_c([1,2,3,4], [1,2,3,4])
    array([ 0.13333333,  0.36666667,  0.66666667,  1.        ,  0.66666667,
            0.36666667,  0.13333333])
    >>> _ncc_c([1,1,1], [1,1,1])
    array([ 0.33333333,  0.66666667,  1.        ,  0.66666667,  0.33333333])
    >>> _ncc_c([1,2,3], [-1,-1,-1])
    array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
    """
    den = np.array(norm(x) * norm(y))
    den[den == 0] = np.Inf

    x_len = len(x)
    fft_size = 1<<(2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
    cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
    return np.real(cc) / den

def reshape(self, vertices):
        """Make snake smoother."""

        arc_length = 0
        # calculate the total arc-length
        for i in range(1, self.length):
            arc_length += norm(vertices[i] - vertices[i - 1])
        # average length for each segment
        seg_length = arc_length / (self.length - 1)
        if self._bc == 'PBC':
            arc_length += norm(vertices[0] - vertices[-1])
            seg_length = arc_length / self.length
        for i in range(1, self.length - 1 if self._bc == 'OBC' else self.length):
            # normalized tangent direction at node i-1
            tan_direction = vertices[i] - vertices[i - 1]
            tan_direction /= norm(tan_direction)
            # move node i
            vertices[i] = vertices[i - 1] + tan_direction * seg_length
        return vertices

def distribution_parameters(self, parameter_name):
        if parameter_name=='trend':
            E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix)
            mean = dot(inv(eye(self.size)+E),self.data)
            cov = (self.parameters.list['sigma2'].current_value)*inv(eye(self.size)+E)
            return {'mean' : mean, 'cov' : cov}
        elif parameter_name=='sigma2':
            E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix)
            pos = self.size
            loc = 0
            scale = 0.5*dot((self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)).T,(self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)))+0.5*dot(dot(self.parameters.list['trend'].current_value.T,E),self.parameters.list['trend'].current_value)
        elif parameter_name=='lambda2':
            pos = self.size-self.total_variation_order-1+self.alpha
            loc = 0.5*(norm(dot(self.derivative_matrix,self.parameters.list['trend'].current_value),ord=1))/self.parameters.list['sigma2'].current_value+self.rho
            scale = 1
        elif parameter_name==str('omega'):
            pos = [sqrt(((self.parameters.list['lambda2'].current_value**2)*self.parameters.list['sigma2'].current_value)/(dj**2)) for dj in dot(self.derivative_matrix,self.parameters.list['trend'].current_value)]
            loc = 0
            scale = self.parameters.list['lambda2'].current_value**2
        return {'pos' : pos, 'loc' : loc, 'scale' : scale}

def initialize_vectors(self):
        self.vec2e = np.zeros(len(self.vocab_e), dtype=np.int)
        for i in range(len(self.vocab_e)):
            w = self.vocab_e[i]
            self.e2vec[w] = i
            self.vec2e[i] = i
        #randomize and normalize the initial vector
        self.vec_e = RD.randn(len(self.vocab_e), self.dim)
        for i in range(len(self.vec_e)):
            self.vec_e[i] /= LA.norm(self.vec_e[i])
        self.vec2r = np.zeros(len(self.vocab_r), dtype=np.int)
        for i in range(len(self.vocab_r)):
            w = self.vocab_r[i]
            #randomize and normalize the initial vector
            self.r2vec[w] = i
            self.vec2r[i] = i
        self.vec_r = RD.randn(len(self.vocab_r), self.dim)
        for i in range(len(self.vec_r)):
            self.vec_r[i] /= LA.norm(self.vec_r[i])
        print "Initialized vectors for all entities and relations."

    #GD of h + r - t

def kNN_relation(self, vec, topk=10, method=0, self_vec_id=None):
        q = []
        for i in range(len(self.vec_r)):
            #skip self
            if self_vec_id != None and i == self_vec_id:
                continue
            if method == 1:
                dist = SP.distance.cosine(vec, self.vec_r[i])
            else:
                dist = LA.norm(vec - self.vec_r[i])
            if len(q) < topk:
                HP.heappush(q, self.index_dist(i, dist))
            else:
                #indeed it fetches the biggest
                tmp = HP.nsmallest(1, q)[0]
                if tmp.dist > dist:
                    HP.heapreplace(q, self.index_dist(i, dist) )
        rst = []
        while len(q) > 0:
            item = HP.heappop(q)
            rst.insert(0, (self.vocab_r[self.vec2r[item.index]], item.dist))
        return rst

    #given relation name, find kNN

def gradient_decent(self,tr,v_e1,v_e2,const_decay, L1=False):
        # f = || v_e1 tr - v_e2 ||_2^2

        assert len(v_e1.shape) == 1
        assert v_e1.shape == v_e2.shape
        assert tr.shape == (v_e1.shape[0], v_e2.shape[0])

        f_res = np.dot(v_e1, tr) - v_e2
        d_v_e1 = 2.0 * np.dot(tr, f_res)
        d_v_e2 = - 2.0 * f_res
        d_tr = 2.0 * np.dot(v_e1[:, np.newaxis], f_res[np.newaxis, :])

        v_e1 -= d_v_e1 * self.rate
        v_e2 -= d_v_e2 * self.rate
        tr -= d_tr * self.rate

        v_e1 /= LA.norm(v_e1)
        v_e2 /= LA.norm(v_e2)
        # don't touch tr

        return LA.norm(np.dot(v_e1, tr) - v_e2)

def kNN_entity_name(self, name, src_lan='en', tgt_lan='fr', topk=10, method=0, replace=True):
        model = self.models.get(src_lan)
        if model == None:
            print "Model for language", src_lan," does not exist."
            return None
        id = model.e2vec.get(name)
        if id == None:
            print name," is not in vocab"
            return None
        if src_lan != tgt_lan:#if you're not quering the kNN in the same language, then no need to get rid of the "self" point. However, transfer the vector.          
            pass_vec = np.dot(model.vec_e[id], self.transfer[src_lan + tgt_lan])
            pass_vec /= LA.norm(pass_vec)
            return self.kNN_entity(pass_vec, tgt_lan, topk, method, self_vec_id=None, replace_q=replace)
        return self.kNN_entity(model.vec_e[id], tgt_lan, topk, method, self_vec_id=id, replace_q=replace)

    #return k nearest relations to a given vector (np.array)

def kNN_relation_name(self, name, src_lan='en', tgt_lan='fr', topk=10, method=0):
        model = self.models.get(src_lan)
        if model == None:
            print "Model for language", src_lan," does not exist."
            return None
        id = model.r2vec.get(name)
        if id == None:
            print name," is not in vocab"
            return None
        if src_lan != tgt_lan:#if you're not quering the kNN in the same language, then no need to get rid of the "self" point    
            pass_vec = np.dot(model.vec_r[id], self.transfer[src_lan + tgt_lan])
            pass_vec /= LA.norm(pass_vec)
            return self.kNN_relation(pass_vec, tgt_lan, topk, method, self_vec_id=None)
        return self.kNN_relation(model.vec_r[id], tgt_lan, topk, method, self_vec_id=id)

    #entity name to vector

def test_matrix_3x3(self):
        # This test has been added because the 2x2 example
        # happened to have equal nuclear norm and induced 1-norm.
        # The 1/10 scaling factor accommodates the absolute tolerance
        # used in assert_almost_equal.
        A = (1 / 10) * \
            np.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
        assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
        assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
        assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
        assert_almost_equal(norm(A, inf), 1.1)
        assert_almost_equal(norm(A, -inf), 0.6)
        assert_almost_equal(norm(A, 1), 1.0)
        assert_almost_equal(norm(A, -1), 0.4)
        assert_almost_equal(norm(A, 2), 0.88722940323461277)
        assert_almost_equal(norm(A, -2), 0.19456584790481812)

def test_bad_args(self):
        # Check that bad arguments raise the appropriate exceptions.

        A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
        B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)

        # Using `axis=<integer>` or passing in a 1-D array implies vector
        # norms are being computed, so also using `ord='fro'`
        # or `ord='nuc'` raises a ValueError.
        assert_raises(ValueError, norm, A, 'fro', 0)
        assert_raises(ValueError, norm, A, 'nuc', 0)
        assert_raises(ValueError, norm, [3, 4], 'fro', None)
        assert_raises(ValueError, norm, [3, 4], 'nuc', None)

        # Similarly, norm should raise an exception when ord is any finite
        # number other than 1, 2, -1 or -2 when computing matrix norms.
        for order in [0, 3]:
            assert_raises(ValueError, norm, A, order, None)
            assert_raises(ValueError, norm, A, order, (0, 1))
            assert_raises(ValueError, norm, B, order, (1, 2))

        # Invalid axis
        assert_raises(ValueError, norm, B, None, 3)
        assert_raises(ValueError, norm, B, None, (2, 3))
        assert_raises(ValueError, norm, B, None, (0, 1, 2))

def signal_by_db(x1, x2, snr, handle_method):
    x1 = x1.astype(np.int32)
    x2 = x2.astype(np.int32)
    l1 = x1.shape[0]
    l2 = x2.shape[0]
    if l1 != l2:
        if handle_method == 'cut':
            ll = min(l1, l2)
            x1 = x1[:ll]
            x2 = x2[:ll]
        elif handle_method == 'append':
            ll = max(l1, l2)
            if l1 < ll:
                x1 = np.append(x1, x1[:ll-l1])
            if l2 < ll:
                x2 = np.append(x2, x2[:ll-l1])

    from numpy.linalg import norm
    x2 = x2 / norm(x2) * norm(x1) / (10.0 ** (0.05 * snr))
    mix = x1 + x2

    return mix

def calc_rotation_matrix(q1, q2, ref_q1, ref_q2):
  ref_nv = np.cross(ref_q1, ref_q2) 
  q_nv = np.cross(q1, q2)
  if min(norm(ref_nv), norm(q_nv)) == 0.:  # avoid 0 degree including angle
    return np.identity(3)
  axis = np.cross(ref_nv, q_nv)
  angle = rad2deg(acos(ref_nv.dot(q_nv) / (norm(ref_nv) * norm(q_nv))))
  R1 = axis_angle_to_rotation_matrix(axis, angle)
  rot_ref_q1, rot_ref_q2 = R1.dot(ref_q1), R1.dot(ref_q2)  # rotate ref_q1,2 plane to q1,2 plane

  cos1 = max(min(q1.dot(rot_ref_q1) / (norm(rot_ref_q1) * norm(q1)), 1.), -1.)  # avoid math domain error
  cos2 = max(min(q2.dot(rot_ref_q2) / (norm(rot_ref_q2) * norm(q2)), 1.), -1.)
  angle1 = rad2deg(acos(cos1))
  angle2 = rad2deg(acos(cos2))
  angle = (angle1 + angle2) / 2.
  axis = np.cross(rot_ref_q1, q1)
  R2 = axis_angle_to_rotation_matrix(axis, angle)

  R = R2.dot(R1)
  return R

def get_lipschitz(data):
    """Get the Lipschitz constant for a specific loss function.

    Only square loss implemented.

    Parameters
    ----------
    data : (n, d) float ndarray
        data matrix
    loss : string
        the selected loss function in {'square', 'logit'}
    Returns
    ----------
    L : float
        the Lipschitz constant
    """
    n, p = data.shape

    if p > n:
        tmp = np.dot(data, data.T)
    else:
        tmp = np.dot(data.T, data)
    return la.norm(tmp, 2)

def prox_l1(w, alpha):
    r"""Proximity operator for l1 norm.

    :math:`\\hat{\\alpha}_{l,m} = sign(u_{l,m})\\left||u_{l,m}| - \\alpha \\right|_+`
    Parameters
    ----------
    u : ndarray
        The vector (of the n-dimensional space) on witch we want
        to compute the proximal operator
    alpha : float
        regularisation parameter
    Returns
    -------
    ndarray : the vector corresponding to the application of the
             proximity operator to u
    """
    return np.sign(w) * np.maximum(np.abs(w) - alpha, 0.)

def __convert_survey_to_sequence(self):
        s = self.__beamline
        if 'LENGTH' not in s:
            s['LENGTH'] = np.nan
        offset = s['ORBIT_LENGTH'][0] / 2.0
        if pd.isnull(offset):
            offset = 0
        self.__beamline['AT_CENTER'] = pd.DataFrame(
            npl.norm(
                [
                    s['X'].diff().fillna(0.0),
                    s['Y'].diff().fillna(0.0)
                ],
                axis=0
            ) - (
                s['LENGTH'].fillna(0.0) / 2.0 - s['ORBIT_LENGTH'].fillna(0.0) / 2.0
            ) + (
                s['LENGTH'].shift(1).fillna(0.0) / 2.0 - s['ORBIT_LENGTH'].shift(1).fillna(0.0) / 2.0
            )).cumsum() / 1000.0 + offset
        self.__converted_from_survey = True

def _initialize(self):
        logger.debug("Initializing Policy.")
        # check if policy is already initialized by the user
        if self.policy.initialized:
            logger.debug("Use pre-set policy parameters.")
            return self.policy.parameters

        # outerwise draw an element at random from the parameter space
        parameter = self.parameter_space.sample()

        for _ in range(1000):
            self.policy.parameters = parameter
            grad = self.estimator(self.policy)

            if (norm(grad) >= 1000 * self.eps):
                return parameter

            parameter = self.parameter_space.sample()

        logger.error('Unable to find non-zero gradient.')

def is_feasible(x, u):
    # Reject going too fast
    for i, v in enumerate(x[3:]):
        if v > velmax_pos_plan[i] or v < velmax_neg_plan[i]:
            return False
    # Body to world
    R = np.array([
                  [np.cos(x[2]), -np.sin(x[2])],
                  [np.sin(x[2]),  np.cos(x[2])],
                ])
    # Boat vertices in world frame
    verts = x[:2] + R.dot(vps).T
    # Check for collisions over all obstacles
    for ob in obs:
        if np.any(npl.norm(verts - ob[:2], axis=1) <= ob[2]):
            return False
    return True

################################################# HEURISTICS

def test_matrix_3x3(self):
        # This test has been added because the 2x2 example
        # happened to have equal nuclear norm and induced 1-norm.
        # The 1/10 scaling factor accommodates the absolute tolerance
        # used in assert_almost_equal.
        A = (1 / 10) * \
            np.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
        assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
        assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
        assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
        assert_almost_equal(norm(A, inf), 1.1)
        assert_almost_equal(norm(A, -inf), 0.6)
        assert_almost_equal(norm(A, 1), 1.0)
        assert_almost_equal(norm(A, -1), 0.4)
        assert_almost_equal(norm(A, 2), 0.88722940323461277)
        assert_almost_equal(norm(A, -2), 0.19456584790481812)

def test_bad_args(self):
        # Check that bad arguments raise the appropriate exceptions.

        A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
        B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)

        # Using `axis=<integer>` or passing in a 1-D array implies vector
        # norms are being computed, so also using `ord='fro'`
        # or `ord='nuc'` raises a ValueError.
        assert_raises(ValueError, norm, A, 'fro', 0)
        assert_raises(ValueError, norm, A, 'nuc', 0)
        assert_raises(ValueError, norm, [3, 4], 'fro', None)
        assert_raises(ValueError, norm, [3, 4], 'nuc', None)

        # Similarly, norm should raise an exception when ord is any finite
        # number other than 1, 2, -1 or -2 when computing matrix norms.
        for order in [0, 3]:
            assert_raises(ValueError, norm, A, order, None)
            assert_raises(ValueError, norm, A, order, (0, 1))
            assert_raises(ValueError, norm, B, order, (1, 2))

        # Invalid axis
        assert_raises(ValueError, norm, B, None, 3)
        assert_raises(ValueError, norm, B, None, (2, 3))
        assert_raises(ValueError, norm, B, None, (0, 1, 2))

def test_polygon(self):
        pts = [[1,0], [1,1], [0,1], [0,2], [6,2]]
        c = CurveFactory.polygon(pts)
        expected_knots = [0,0,1,2,3,9,9]
        actual_knots   = c.knots(0,True)

        self.assertEqual(len(c),      5)
        self.assertEqual(c.order(0),  2)
        self.assertEqual(c.dimension, 2)
        self.assertAlmostEqual(np.linalg.norm(expected_knots - actual_knots), 0.0)

        c = CurveFactory.polygon([0,0], [1,0], [0,1], [-1,0], relative=True)
        self.assertEqual(len(c), 4)
        self.assertAlmostEqual(c[2][0], 1)
        self.assertAlmostEqual(c[2][1], 1)
        self.assertAlmostEqual(c[3][0], 0)
        self.assertAlmostEqual(c[3][1], 1)

def test_bezier(self):
        crv = CurveFactory.bezier([[0,0], [0,1], [1,1], [1,0], [2,0], [2,1],[1,1]])
        self.assertEqual(len(crv.knots(0)), 3)
        self.assertTrue(np.allclose(crv(0), [0,0]))
        t = crv.tangent(0)
        self.assertTrue(np.allclose(t/norm(t), [0,1]))
        t = crv.tangent(1, above=False)
        self.assertTrue(np.allclose(t/norm(t), [0,-1]))
        t = crv.tangent(1, above=True)
        self.assertTrue(np.allclose(t/norm(t), [1,0]))
        self.assertTrue(np.allclose(crv(1), [1,0]))
        self.assertTrue(crv.order(0), 4)

        # test the exact same curve, only with relative keyword
        crv = CurveFactory.bezier([[0,0], [0,1], [1,0], [0,-1], [1,0], [0,1],[1,0]], relative=True)
        self.assertEqual(len(crv.knots(0)), 3)
        self.assertTrue(np.allclose(crv(0), [0,0]))
        t = crv.tangent(0)
        self.assertTrue(np.allclose(t/norm(t), [0,1]))
        t = crv.tangent(1, above=False)
        self.assertTrue(np.allclose(t/norm(t), [0,-1]))
        t = crv.tangent(1, above=True)
        self.assertTrue(np.allclose(t/norm(t), [1,0]))
        self.assertTrue(np.allclose(crv(1), [1,0]))
        self.assertTrue(crv.order(0), 4)

def test_volume_sphere(self):
        # unit ball
        ball = VolumeFactory.sphere(type='square')
        # test 7x7 grid for radius = 1
        for surf in ball.faces():
            for u in np.linspace(surf.start(0), surf.end(0), 7):
                for v in np.linspace(surf.start(1), surf.end(1), 7):
                    self.assertAlmostEqual(np.linalg.norm(surf(u, v), 2), 1.0)
            self.assertAlmostEqual(surf.area(), 4*pi/6, places=3)

        # unit ball
        ball = VolumeFactory.sphere(type='radial')
        # test 7x7 grid for radius = 1
        for u in np.linspace(surf.start(0), surf.end(0), 7):
            for v in np.linspace(surf.start(1), surf.end(1), 7):
                self.assertAlmostEqual(np.linalg.norm(ball(u, v, 0), 2), 1.0) # w=0 is outer shell
                self.assertAlmostEqual(np.linalg.norm(ball(u, v, 1), 2), 0.0) # w=1 is the degenerate core
        self.assertAlmostEqual(ball.faces()[4].area(), 4*pi, places=3)

def test_cylinder_surface(self):
        # unit cylinder
        surf = SurfaceFactory.cylinder()
        # test 7x7 grid for xy-radius = 1 and v=z
        for u in np.linspace(surf.start()[0], surf.end()[0], 7):
            for v in np.linspace(surf.start()[1], surf.end()[1], 7):
                x = surf(u, v)
                self.assertAlmostEqual(np.linalg.norm(x[:2], 2), 1.0)  # (x,y) coordinates to z-axis
                self.assertAlmostEqual(x[2], v)  # z coordinate should be linear
        self.assertAlmostEqual(surf.area(), 2*pi, places=3)

        # cylinder with parameters
        surf = SurfaceFactory.cylinder(r=2, h=5, center=(0,0,1), axis=(1,0,0))
        for u in np.linspace(surf.start()[0], surf.end()[0], 7):
            for v in np.linspace(surf.start()[1], surf.end()[1], 7):
                x = surf(u, v)
                self.assertAlmostEqual(x[1]**2+(x[2]-1)**2, 2.0**2) # distance to (z-1)=y=0
                self.assertAlmostEqual(x[0], v*5)                   # x coordinate should be linear
        self.assertAlmostEqual(surf.area(), 2*2*pi*5, places=3)

def test_cylinder_volume(self):
        # unit cylinder
        vol = VolumeFactory.cylinder()
        # test 5x5x5 grid for xy-radius = w and v=z
        for u in np.linspace(vol.start()[0], vol.end()[0], 5):
            for v in np.linspace(vol.start()[1], vol.end()[1], 5):
                for w in np.linspace(vol.start()[2], vol.end()[2], 5):
                    x = vol(u, v, w)
                    self.assertAlmostEqual(
                        np.linalg.norm(x[:2], 2), u)  # (x,y) coordinates to z-axis
                    self.assertAlmostEqual(x[2], w)  # z coordinate should be linear
        self.assertAlmostEqual(vol.volume(), pi, places=3)

        # cylinder with parameters
        vol = VolumeFactory.cylinder(r=2, h=5, center=(0,0,1), axis=(1,0,0))
        for u in np.linspace(vol.start()[0], vol.end()[0], 5):
            for v in np.linspace(vol.start()[1], vol.end()[1], 5):
                for w in np.linspace(vol.start()[2], vol.end()[2], 5):
                    x = vol(u, v, w)
                    self.assertAlmostEqual(x[1]**2+(x[2]-1)**2, u**2)
                    self.assertAlmostEqual(x[0], w*5)
        self.assertAlmostEqual(vol.volume(), 2*2*pi*5, places=3)

def eval(self, coords, grad=False):
        v1 = (coords[self.i]-coords[self.j])/bohr
        v2 = (coords[self.k]-coords[self.j])/bohr
        dot_product = np.dot(v1, v2)/(norm(v1)*norm(v2))
        if dot_product < -1:
            dot_product = -1
        elif dot_product > 1:
            dot_product = 1
        phi = np.arccos(dot_product)
        if not grad:
            return phi
        if abs(phi) > pi-1e-6:
            grad = [
                (pi-phi)/(2*norm(v1)**2)*v1,
                (1/norm(v1)-1/norm(v2))*(pi-phi)/(2*norm(v1))*v1,
                (pi-phi)/(2*norm(v2)**2)*v2
            ]
        else:
            grad = [
                1/np.tan(phi)*v1/norm(v1)**2-v2/(norm(v1)*norm(v2)*np.sin(phi)),
                (v1+v2)/(norm(v1)*norm(v2)*np.sin(phi)) -
                1/np.tan(phi)*(v1/norm(v1)**2+v2/norm(v2)**2),
                1/np.tan(phi)*v2/norm(v2)**2-v1/(norm(v1)*norm(v2)*np.sin(phi))
            ]
        return phi, grad

def quadratic_step(g, H, w, trust, log=no_log):
    ev = np.linalg.eigvalsh((H+H.T)/2)
    rfo = np.vstack((np.hstack((H, g[:, None])),
                     np.hstack((g, 0))[None, :]))
    D, V = np.linalg.eigh((rfo+rfo.T)/2)
    dq = V[:-1, 0]/V[-1, 0]
    l = D[0]
    if norm(dq) <= trust:
        log('Pure RFO step was performed:')
        on_sphere = False
    else:
        def steplength(l):
            return norm(np.linalg.solve(l*eye(H.shape[0])-H, g))-trust
        l = Math.findroot(steplength, ev[0])  # minimization on sphere
        dq = np.linalg.solve(l*eye(H.shape[0])-H, g)
        on_sphere = False
        log('Minimization on sphere was performed:')
    dE = dot(g, dq)+0.5*dq.dot(H).dot(dq)  # predicted energy change
    log('* Trust radius: {:.2}'.format(trust))
    log('* Number of negative eigenvalues: {}'.format((ev < 0).sum()))
    log('* Lowest eigenvalue: {:.3}'.format(ev[0]))
    log('* lambda: {:.3}'.format(l))
    log('Quadratic step: RMS: {:.3}, max: {:.3}'.format(Math.rms(dq), max(abs(dq))))
    log('* Predicted energy change: {:.3}'.format(dE))
    return dq, dE, on_sphere

def XB(Y,index_PQ, index_P, n_PQ, n_P):
    case_number, _ = np.shape(Y)
    G, B = iterator.util.get_G_B(Y)
    X = np.zeros((case_number, case_number))
    B_p = np.zeros((n_P,n_P))
    B_pp = np.zeros((n_PQ,n_PQ))
    for i in xrange(case_number):
        for j in xrange(case_number):
            if LA.norm(Y[i][j]) > 1e-5 and i != j:
                X[i][j] = np.reciprocal(np.imag(np.reciprocal(Y[i][j])))

    for i in xrange(case_number):
        for j in xrange(case_number):
            if i != j:
                X[i][i] -= X[i][j]

    for i in xrange(0, n_P):
        for j in xrange(0, n_P):
            B_p[i][j] = X[index_P[i]][index_P[j]]
    #---------------------------------------------------
    for i in xrange(0, n_PQ):
        for j in xrange(0, n_PQ):
            B_pp[i][j] = B[index_PQ[i]][index_PQ[j]]

    return B_p, B_pp

def predict(self, X):
        X = self._pre_processing_x(X)
        Y = np.zeros((X.shape[0], self.n_class))
        A = self.A

        # because X is (N x p), A is (K x p), we can to get the X_star (NxK)
        X_star = X @ A.T

        for k in range(self.n_class):
            # mu_s_star shape is (p,)
            mu_k_star = A @ self.Mu[k]

            # Ref: http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html
            # Ref: http://stackoverflow.com/questions/1401712/how-can-the-euclidean-distance-be-calculated-with-numpy
            Y[:, k] = LA.norm(X_star - mu_k_star, axis=1) * 0.5 - log(self.Pi[k])

        # Python index start from 0, transform to start with 1
        y_hat = Y.argmin(axis=1).reshape((-1, 1)) + 1
        return y_hat

def normalise_word_vectors(word_vectors, norm=1.0):
    """
    This method normalises the collection of word vectors provided in the word_vectors dictionary.
    """
    for word in word_vectors:
        word_vectors[word] /= math.sqrt((word_vectors[word]**2).sum() + 1e-6)
        word_vectors[word] = word_vectors[word] * norm
    return word_vectors

def normalise_vector(v1):
    return v1 / norm(v1)

def distance(v1, v2, normalised_vectors=False):
    """
    Returns the cosine distance between two vectors. 
    If the vectors are normalised, there is no need for the denominator, which is always one. 
    """
    if normalised_vectors:
        return 1 - dot(v1, v2)
    else:
        return 1 - dot(v1, v2) / ( norm(v1) * norm(v2) )

def Verify(**kwargs):
    msg, A, m, n, sd, q, eta, z, c, kappa = kwargs['msg'], kwargs['A'], kwargs['m'], kwargs['n'], kwargs['sd'], kwargs['q'], kwargs['eta'], kwargs['z'], kwargs['c'], kwargs['kappa']
    B2 = eta*sd*np.sqrt(m)
    reduced_prod = util.vector_to_Zq(np.matmul(A,z) + q*c, 2*q)
    #print np.sqrt(z.dot(z)),B2
    #print LA.norm(z,np.inf),float(q)/4
    if np.sqrt(z.dot(z)) > B2  or LA.norm(z,np.inf) >= float(q)/4:      
        return False    
    if np.array_equal(c, hash_iterative(np.array_str(reduced_prod)+msg, n, kappa)):
        return True
    return False

def normal_direction(p, p1, p2):
        """Compute normal direction at p in the segment p1->p->p2."""

        e1 = (p1 - p) / norm(p1 - p)
        e2 = Snake.rotate(e1, np.pi / 2.)
        x = np.dot(p2 - p, e1)
        y = np.dot(p2 - p, e2)
        theta = np.arctan2(y, x)
        return Snake.rotate(e1, theta / 2.)

def vector_correlation(vect1,vect2):
    """
    Compute correlation between two vector, which is the the cosine of the angle between two vectors in Euclidean space of any number of dimensions.
    The dot product is directly related to the cosine of the angle between two vectors if they are normed !!!
    """
    # -- We make sure that we have normed vectors.
    from numpy.linalg import norm
    if (np.round(norm(vect1)) != 1.):
        vect1 = vect1/norm(vect1)
    if (np.round(norm(vect2)) != 1.):
        vect2 = vect2/norm(vect2)

    return np.round(np.dot(vect1,vect2),3)

def admm_phase2(x0, prob, rho, tol=1e-2, num_iters=1000, viol_lim=1e4):
    logging.info("Starting ADMM phase 2 with rho %.3f", rho)

    bestx = np.copy(x0)

    z = np.copy(x0)
    xs = [np.copy(x0) for i in range(prob.m)]
    us = [np.zeros(prob.n) for i in range(prob.m)]

    if prob.rho != rho:
        prob.rho = rho
        zlhs = 2*(prob.f0.P + rho*prob.m*sp.identity(prob.n)).tocsc()
        prob.z_solver = SLA.factorized(zlhs)

    last_z = None
    for t in range(num_iters):
        rhs = 2*rho*(sum(xs)-sum(us)) - prob.f0.qarray
        z = prob.z_solver(rhs)

        # TODO: parallel x/u-updates
        for i in range(prob.m):
            xs[i] = onecons_qcqp(z + us[i], prob.fi(i))
        for i in range(prob.m):
            us[i] += z - xs[i]

        # TODO: termination condition
        if last_z is not None and LA.norm(last_z - z) < tol:
            break
        last_z = z

        maxviol = max(prob.violations(z))
        logging.info("Iteration %d, violation %.3f", t, maxviol)

        if maxviol > viol_lim: break
        bestx = np.copy(prob.better(z, bestx))

    return bestx

def vector_norm(data, axis=None, out=None):
    """Return length, i.e. Euclidean norm, of ndarray along axis.

    >>> v = numpy.random.random(3)
    >>> n = vector_norm(v)
    >>> numpy.allclose(n, numpy.linalg.norm(v))
    True
    >>> v = numpy.random.rand(6, 5, 3)
    >>> n = vector_norm(v, axis=-1)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=2)))
    True
    >>> n = vector_norm(v, axis=1)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
    True
    >>> v = numpy.random.rand(5, 4, 3)
    >>> n = numpy.empty((5, 3))
    >>> vector_norm(v, axis=1, out=n)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
    True
    >>> vector_norm([])
    0.0
    >>> vector_norm([1])
    1.0

    """
    data = numpy.array(data, dtype=numpy.float64, copy=True)
    if out is None:
        if data.ndim == 1:
            return math.sqrt(numpy.dot(data, data))
        data *= data
        out = numpy.atleast_1d(numpy.sum(data, axis=axis))
        numpy.sqrt(out, out)
        return out
    else:
        data *= data
        numpy.sum(data, axis=axis, out=out)
        numpy.sqrt(out, out)

def cos_distance(cls, x, y ):
        x, y = np.mat(x), np.mat(y)
        #num = float(x.dot(y.T))
        num = float(x * y.T)
        denom = la.norm(x) * la.norm(y)
        #dist = 0.5 + 0.5 * (num / denom)
        dist = 1.0 if denom == 0 else num / denom 
        return 3 * (1 - dist)

def save_problem(base,prob):
    print('saving {b}.mat,{b}.npz norm(x)={x:.7f} norm(y)={y:.7f}'.format(b=base,x=la.norm(prob.xval), y=la.norm(prob.yval) ) )
    D=dict(A=prob.A,x=prob.xval,y=prob.yval)
    np.savez( base + '.npz', **D)
    savemat(base + '.mat',D,oned_as='column')

def bernoulli_gaussian_trial(M=250,N=500,L=1000,pnz=.1,kappa=None,SNR=40):

    A = np.random.normal(size=(M, N), scale=1.0 / math.sqrt(M)).astype(np.float32)
    if kappa >= 1:
        # create a random operator with a specific condition number
        U,_,V = la.svd(A,full_matrices=False)
        s = np.logspace( 0, np.log10( 1/kappa),M)
        A = np.dot( U*(s*np.sqrt(N)/la.norm(s)),V).astype(np.float32)
    A_ = tf.constant(A,name='A')
    prob = TFGenerator(A=A,A_=A_,pnz=pnz,kappa=kappa,SNR=SNR)
    prob.name = 'Bernoulli-Gaussian, random A'

    bernoulli_ = tf.to_float( tf.random_uniform( (N,L) ) < pnz)
    xgen_ = bernoulli_ * tf.random_normal( (N,L) )
    noise_var = pnz*N/M * math.pow(10., -SNR / 10.)
    ygen_ = tf.matmul( A_,xgen_) + tf.random_normal( (M,L),stddev=math.sqrt( noise_var ) )

    prob.xval = ((np.random.uniform( 0,1,(N,L))<pnz) * np.random.normal(0,1,(N,L))).astype(np.float32)
    prob.yval = np.matmul(A,prob.xval) + np.random.normal(0,math.sqrt( noise_var ),(M,L))
    prob.xinit = ((np.random.uniform( 0,1,(N,L))<pnz) * np.random.normal(0,1,(N,L))).astype(np.float32)
    prob.yinit = np.matmul(A,prob.xinit) + np.random.normal(0,math.sqrt( noise_var ),(M,L))
    prob.xgen_ = xgen_
    prob.ygen_ = ygen_
    prob.noise_var = noise_var

    return prob

def add_noise(self,Y0):
        'add noise at the given SNR, returns Y0+W,wvar'
        wvar = (la.norm(Y0)**2/Y0.size) * 10**(-self.SNR_dB/10)
        if self.cpx:
            Y =(Y0 + crandn(Y0.shape) * sqrt(wvar/2)).astype(np.complex64,copy=False)
        else:
            Y = (Y0 + np.random.normal(scale=sqrt(wvar),size=Y0.shape) ).astype(np.float32,copy=False)
        return Y,wvar

def _test_awgn(self,cpx):
            snr = np.random.uniform(3,20)
            p = Problem(cpx=cpx,SNR_dB=snr)
            X = p.genX(5)
            self.assertEqual( np.iscomplexobj(X) , cpx )
            Y0 = p.fwd(X)
            self.assertEqual( np.iscomplexobj(Y0) , cpx )
            Y,wvar = p.add_noise(Y0)
            self.assertEqual( np.iscomplexobj(Y) , cpx )
            snr_obs = -20*np.log10( la.norm(Y-Y0)/la.norm(Y0))
            self.assertTrue( abs(snr-snr_obs) < 1.0, 'gross error in add_noise')
            wvar_obs = la.norm(Y0-Y)**2/Y.size
            self.assertTrue( .5 < wvar_obs/wvar < 1.5, 'gross error in add_noise wvar')