def _rank1(X, profunc1, args1, profunc2, args2, niter=500, eps=1e-6):
	u = X.sum(axis=1); u /= norm(u)
	v = X.sum(axis=0); v /= norm(v)

	rho_old = dot(dot(u,X),v)
	for i in xrange(niter):
		alpha = dot(X,v)
		if profunc1 == None:
			u = alpha
		else:
			u = profunc1(alpha, **args1)
		u /= norm(u)

		beta = dot(X.T,u)
		if profunc2 == None:
			v = beta
		else:
			v = profunc2(beta, **args2)
		rho = norm(v)
		v /= rho

		if abs(rho-rho_old) <= eps:
			break
		else:
			rho_old = rho

	return u, rho, v 

def isparallel(O1,O2):
    '''
    Judge whether two array-like vectors are parallel to each other.

    Parameters
    ----------
    O1,O2 : 1d array-like
        The input vectors.

    Returns
    -------
    int
        *  0: not parallel
        *  1: parallel
        * -1: anti-parallel
    '''
    norm1=nl.norm(O1)
    norm2=nl.norm(O2)
    if norm1<RZERO or norm2<RZERO:
        return 1
    elif O1.shape[0]==O2.shape[0]:
        buff=np.inner(O1,O2)/(norm1*norm2)
        if np.abs(buff-1)<RZERO:
            return 1
        elif np.abs(buff+1)<RZERO:
            return -1
        else:
            return 0
    else:
        raise ValueError("isparallel error: the shape of the array-like vectors does not match.")

def test3():
    import spacy.en
    from spacy.parts_of_speech import ADV
    nlp = spacy.en.English()
    # Find log probability of Nth most frequent word
    probs = [lex.prob for lex in nlp.vocab]
    probs.sort()
    is_adverb = lambda tok: tok.pos == ADV and tok.prob < probs[-1000]
    tokens = nlp(u"‘Give it back,’ he pleaded abjectly, ‘it’s mine.’")
    o = u''.join(tok.string.upper() if is_adverb(tok) else tok.string for tok in tokens)
    assert o == u'‘Give it back,’ he pleaded ABJECTLY, ‘it’s mine.’'

    pleaded = tokens[7]
    assert pleaded.repvec.shape == (300,)
    o = pleaded.repvec[:5]
    assert sum(o) != 0
    from numpy import dot
    from numpy.linalg import norm

    cosine = lambda v1, v2: dot(v1, v2) / (norm(v1) * norm(v2))
    words = [w for w in nlp.vocab if w.check(IS_LOWER) and w.has_repvec]
    words.sort(key=lambda w: cosine(w.repvec, pleaded.repvec))
    words.reverse()
    o = [w.orth_ for w in words[0:20]]
    assert o == [u'pleaded', u'pled', u'plead', u'confessed', u'interceded',
                 u'pleads', u'testified', u'conspired', u'motioned', u'demurred',
                 u'countersued', u'remonstrated', u'begged', u'apologised',
                 u'consented', u'acquiesced', u'petitioned', u'quarreled',
                 u'appealed', u'pleading']
    o = [w.orth_ for w in words[50:60]]
    assert o == [u'counselled', u'bragged', u'backtracked', u'caucused', u'refiled',
                 u'dueled', u'mused', u'dissented', u'yearned', u'confesses']
    o = [w.orth_ for w in words[100:110]]
    assert o == [u'cabled', u'ducked', u'sentenced', u'perjured', u'absconded',
                 u'bargained', u'overstayed', u'clerked', u'confided', u'sympathizes']

def sparse_calculate_error(mat, sketch, normalized=True):
    cov =  (mat.transpose()).dot(mat).todense()
    cov_sketch = np.dot(sketch.T, sketch)
    if normalized:
        return (ln.norm(cov - cov_sketch, ord=2) / sparse_squared_frobenius_norm(mat))
    else:
        return ln.norm(cov - cov_sketch, ord=2)

    def test_retrieve_data(self):
        ptree = PropertyTree()
        ptree.put_string('type', 'SeriesRC')
        ptree.put_double('series_resistance', 100e-3)
        ptree.put_double('capacitance', 2.5)
        device = EnergyStorageDevice(ptree)

        ptree = PropertyTree()
        ptree.put_string('type', 'ElectrochemicalImpedanceSpectroscopy')
        ptree.put_double('frequency_upper_limit', 1e+2)
        ptree.put_double('frequency_lower_limit', 1e-1)
        ptree.put_int('steps_per_decade', 1)
        ptree.put_int('steps_per_cycle', 64)
        ptree.put_int('cycles', 2)
        ptree.put_int('ignore_cycles', 1)
        ptree.put_double('dc_voltage', 0)
        ptree.put_string('harmonics', '3')
        ptree.put_string('amplitudes', '5e-3')
        ptree.put_string('phases', '0')
        eis = Experiment(ptree)

        with File('trash.hdf5', 'w') as fout:
            eis.run(device, fout)
        spectrum_data = eis._data

        with File('trash.hdf5', 'r') as fin:
            retrieved_data = retrieve_impedance_spectrum(fin)

        print(spectrum_data['impedance'] - retrieved_data['impedance'])
        print(retrieved_data)
        self.assertEqual(linalg.norm(spectrum_data['frequency'] -
                                     retrieved_data['frequency'], inf), 0.0)
        # not sure why we don't get equality for the impedance
        self.assertLess(linalg.norm(spectrum_data['impedance'] -
                                    retrieved_data['impedance'], inf), 1e-10)

    def evaluate(self, locations):
        total = 0
        for i, layer in enumerate(self._all_edges):
            s = layer.shape
            for j1 in range(s[0]):
                if j1 == self.nodes[i].shape[1]:
                    break
                for j2 in range(s[1]):
                    #firstpoint = (i*hscale + hoff, j1*vscale + voff)
                    #secondpoint = ((i+1)*hscale + hoff, j2*vscale + voff)
                    one = locations[(i, j1)]
                    two = locations[(i + 1, j2)]
                    distance = linalg.norm(one - two)
                    f = distance * distance

                    influence = abs(layer[j1, j2])
                    # quality is decreased when connected things are far apart
                    total -= f * influence

        for _, v in locations:
            for _, v2 in locations:
                if v != v2:
                    d = linalg.norm(v - v2)
                    total += 1000 / d  # this is like electric potential

        return total

def compute_TPS_K(ctrl_pts, landmarks = None, _lambda = 0):
    """
    Compute the kernel matrix for thin-plate splines.
    Reference:
      Landmark-based Image Analysis, Karl Rohr, p195
    """

    #kernel_func = [lambda r,_lambda=0: 0 if r==0 else r*r*log(r), lambda r,_lambda=0: -r]
    #the above definition is not used because the if else syntax is not supported in Python2.4
    def kernel_func_2d(r, _lambda=0):
        #_lambda reserved for regularization
        if r == 0:
            return 0
        else:
            return r*r*log(r)
    def kernel_func_3d(r, _lambda=0):
        #_lambda reserved for regularization
        return -r
    kernel_func = (kernel_func_2d, kernel_func_3d)

    [n,d] = ctrl_pts.shape
    K = [kernel_func[d-2](norm(ctrl_pts[i]-ctrl_pts[j]), _lambda) for i in arange(n) for j in arange(n)]
    K = array(K).reshape(n,n)
    if landmarks is not None:
        [m,d] = landmarks.shape  # assert (d,d) equal
        U = [kernel_func[d-2](norm(landmarks[i]-ctrl_pts[j]), _lambda) for i in arange(m) for j in arange(n)]
        U = array(U).reshape(m,n)
    else:
        U = None
    return K,U

    def distance(self, other, method='euclidean'):
        """
        Distance between the center of this source and another.

        Parameters
        ----------
        other : Source, or array-like
            Either another source, or the center coordinates of another source

        method : str
            Specify a distance measure to used for spatial distance between source
            centers. Current options include Euclidean distance ('euclidean') and 
            L1-norm ('l1'). 

        """
        from numpy.linalg import norm

        checkParams(method, ['euclidean', 'l1'])

        if method == 'l1':
            order = 1
        else:
            order = 2

        if isinstance(other, Source):
            return norm(self.center - other.center, ord=order)
        elif isinstance(other, list) or isinstance(other, ndarray):
            return norm(self.center - asarray(other), ord=order)

    def reset(self, direction, speed):
        # initial location
        p = [random.uniform(-1, 1) for i in range(3)]
        p /= norm(p)

        # orienting point
        o = [0, 0, 0]
        mini = min(range(len(p)), key=lambda i: abs(p[i]))
        o[mini] = 1 if p[mini] < 0 else -1

        # velocity vector
        v = cross(p, o)
        v /= norm(v)
        v *= speed*pi/180

        r = 1.145
        shape = SphericalPolygon([rotate(rotate(p, o, r*random.uniform(0.9,1.1)), p, th)
                                  for th in [i*pi/8 for i in range(16)]])

        for t in self.tiles.values():
            t.bottom = 0
            t.layers = [Layer('T', 1)] if shape.contains(t.vector) else []
            t.limit()

        self._indexedtiles = []
        for t in self.tiles.values():
            self._indexedtiles.append(t)

        self._index = PointTree(dict([[self._indexedtiles[i].vector, i]
                                      for i in range(len(self._indexedtiles))]))
    
        self._direction = direction
        self._velocity = v

 def test_ordinary(self):
     N=4
     np.random.seed()
     state,target=np.zeros((2,)*N),SQN(0.0)
     for index in QuantumNumbers.decomposition([SQNS(0.5)]*N,signs=[1]*N,target=target):
         state[index]=np.random.random()
     state=state.reshape((-1,))
     sites=[Label('S%s'%i,qns=SQNS(0.5),flow=1) for i in range(N)]
     bonds=[Label('B%s'%i,qns=None,flow=None) for i in range(N+1)]
     bonds[+0]=bonds[+0].replace(qns=SQNS(0.0),flow=+1)
     bonds[-1]=bonds[-1].replace(qns=QuantumNumbers.mono(target),flow=-1)
     for cut in range(N+1):
         mps=MPS.fromstate(state,sites,bonds,cut=cut,ttype='D')
         self.assertTrue(all(mps.iscanonical()))
         self.assertAlmostEqual(norm(state-mps.state),0.0)
     for cut in range(N+1):
         mps.canonicalize(cut)
         self.assertTrue(all(mps.iscanonical()))
     for cut in range(N+1):
         mps=MPS.fromstate(state,sites,bonds,cut=cut,ttype='S')
         self.assertTrue(all(mps.iscanonical()))
         self.assertAlmostEqual(norm(state-mps.state),0.0)
     for cut in range(N+1):
         mps.canonicalize(cut)
         self.assertTrue(all(mps.iscanonical()))

def cosine_similarity(a, b):
    tn = np.inner(a, b)
    td = la.norm(a) * la.norm(b)
    if td != 0.0:
        return tn / td
    else:
        return 0.0

def summarize_evaluation(query=None, url=None, summary=None):
    j=[]
    if url:
        b = URL(url)
        a = Document(b.download(cached=True))
        for b in a.get_elements_by_tagname("p"):
            j.append(plaintext(b.content).encode("utf-8"))
        j = [word for sentence in j for word in sentence.split() if re.match("^[a-zA-Z_-]*$", word) or '.' in word or "'" in word or '"' in word]
        j = ' '.join(j)
        lsa = LSA(stopwords, ignore_characters)
        sentences = j.split('.')
        sentences = [sentence for sentence in sentences if len(sentence)>1 and sentence != '']
        for sentence in sentences:
            lsa.parse(sentence)
    else:
        lsa = LSA(stopwords, ignore_characters)
        for sentence in query:
            lsa.parse(sentence)
    lsa.build()
    lsa.calc()
    lsa2 = LSA(stopwords, ignore_characters)
    for sentence in summary:
        lsa2.parse(sentence)
    lsa2.build()
    lsa2.calc()
    vectors =[(dot(lsa.S,lsa.U[0,:]),dot(lsa.S,lsa.U[i,:])) for i in range(len(lsa.U))]
    vectors2 =[(dot(lsa2.S,lsa2.U[0,:]),dot(lsa2.S,lsa2.U[i,:])) for i in range(len(lsa2.U))]
    angles = [arccos(dot(a,b)/(norm(a,2)*norm(b,2))) for a in vectors for b in vectors2]
    return str(abs(1 - float(angles[1])/float(pi/2)))

def check_cost_function(lambda_coef=0):
    X_t = np.random.rand(4, 3)
    Theta_t = np.random.rand(5, 3)
    Y = X_t.dot(Theta_t.T)
    Y[np.random.rand(*Y.shape) > 0.5] = 0

    R = np.zeros_like(Y)
    R[Y != 0] = 1

    X = np.random.randn(*X_t.shape);
    Theta = np.random.randn(*Theta_t.shape);
    num_movies, num_users = Y.shape
    num_features = Theta_t.shape[1]

    J = lambda t: cofi_cost_func(t, Y, R, num_users, num_movies, num_features, lambda_coef)[0]
    numgrad = compute_numerical_gradient(J, np.hstack((X.ravel(), Theta.ravel())))

    cost, grad = cofi_cost_func(np.hstack((X.ravel(), Theta.ravel())),  Y, R, num_users, num_movies, num_features, lambda_coef)
    for i, j in zip(numgrad, grad):
        print i, j
    print('The above two columns you get should be very similar.\n'
         '(Left-Your Numerical Gradient, Right-Analytical Gradient)\n')
    diff = norm(numgrad-grad)/norm(numgrad+grad)
    print('If your backpropagation implementation is correct, then \n'
         'the relative difference will be small (less than 1e-9). \n'
         'Relative Difference: %s' % diff)

def AreaNormal(nodes):
    """
    Returns area,unitNormal
    n = Normal = a x b
    Area   = 1/2 * |a x b|
    V = <v1,v2,v3>
    |V| = sqrt(v1^0.5+v2^0.5+v3^0.5) = norm(V)

    Area = 0.5 * |n|
    unitNormal = n/|n|
    """
    (n0, n1, n2) = nodes
    a = n0 - n1
    b = n0 - n2
    vector = cross(a, b)
    length = norm(vector)
    normal = vector / length
    area = 0.5 * length
    if not allclose(norm(normal), 1.):
        print("a = ", a)
        print("b = ", b)
        print("normal = ", normal)
        print("length = ", length)
        sys.exit('check...')
    return area, normal

def hexagon_bz(reciprocals=None,nk=100,vh='H'):
    '''
    The whole hexagonal BZ.
    '''
    if reciprocals is not None:
        b1=reciprocals[0]
        b2=reciprocals[1]
        temp=np.inner(b1,b2)/nl.norm(b1)/nl.norm(b2)
        assert np.abs(np.abs(temp)-0.5)<RZERO
        if np.abs(temp+0.5)<RZERO: b2=-b2
    else:
        if vh in ('H','h'):
            b1=np.array([np.sqrt(3.0)/2,0.5])*4*np.pi/np.sqrt(3.0)
            b2=np.array([np.sqrt(3.0)/2,-0.5])*4*np.pi/np.sqrt(3.0)
        else:
            b1=np.array([1.0,0.0])*4*np.pi/np.sqrt(3.0)
            b2=np.array([0.5,np.sqrt(3.0)/2])*4*np.pi/np.sqrt(3.0)
    p0,p1,p2,p3,p4=-(b1+b2)/3,(b1+b2)/3,(b1+b2)*2/3,(b1*2-b2)/3,(b2*2-b1)/3
    mesh=np.zeros((nk**2,b1.shape[0]))
    for i in range(nk):
        for j in range(nk):
            coords=b1*i/nk+b2*j/nk+p0
            if isintratriangle(coords,p1,p2,p3,vertexes=(False,True,False),edges=(True,True,False)): coords=coords-b1
            if isintratriangle(coords,p1,p2,p4,vertexes=(False,True,False),edges=(True,True,False)): coords=coords-b2
            mesh[i*nk+j,:]=coords
    volume=np.abs(np.cross(b1,b2))
    return BaseSpace(('k',mesh,volume))

def Normal(a, b):
    """finds the unit normal vector of 2 vectors"""
    vector = cross(a, b)
    length = norm(vector)
    normal = vector / length
    assert allclose(norm(normal), 1.)
    return normal

    def __init__(self,reciprocals,path):
        '''
        Constructor.

        Parameters
        ----------
        reciprocals : iterable of 1d ndarray
            The translation vectors of the reciprocal lattice.
        path : str
            The str-formed path.
        '''
        path=path.replace(' ', '')
        assert path[0] in KMap.database and path[1]==':'
        space,path,database,reciprocals=path[0],path[2:].split(','),KMap.database[path[0]],np.asarray(reciprocals)
        if space=='L':
            assert len(reciprocals)==1
        elif space=='S':
            assert len(reciprocals)==2
            inner=np.inner(reciprocals[0],reciprocals[1])/nl.norm(reciprocals[0])/nl.norm(reciprocals[1])
            assert np.abs(inner)<RZERO
        elif space=='H':
            assert len(reciprocals)==2
            inner=np.inner(reciprocals[0],reciprocals[1])/nl.norm(reciprocals[0])/nl.norm(reciprocals[1])
            assert np.abs(np.abs(inner)-0.5)<RZERO
            if np.abs(inner+0.5)<RZERO: reciprocals[1]=-reciprocals[1]
        for segment in path:
            segment=segment.split('-')
            assert len(segment)==2
            self.append([reciprocals.T.dot(database[segment[0]]),reciprocals.T.dot(database[segment[1]])])

def genGraph(S_actual, S_est, S_previous, empCov_set, nodeID, e1, e2, e3, e4, display = False):
    D = np.where(S_est != 0)[0].shape[0]
    T = np.where(S_actual != 0)[0].shape[0]
    TandD = float(np.where(np.logical_and(S_actual,S_est) == True)[0].shape[0])
    P = TandD/D
    R = TandD/T
    offDiagDiff = S_actual - S_est
    offDiagDiff = offDiagDiff - np.diag(np.diag(offDiagDiff))
    S_diff = (S_est - S_previous)  
    S_diff = S_diff - np.diag(np.diag(S_diff))
    ind = (S_diff < 1e-2) & (S_diff > - 1e-2)
    S_diff[ind] = 0    
    K = np.count_nonzero(S_diff)
    e1.append( alg.norm(offDiagDiff, 'fro'))
    e2.append(2* P*R/(P+R))
    
    
    K = float(np.where(np.logical_and((S_est>0) != (S_previous>0), S_est>0) == True)[0].shape[0])
    e3.append(-np.log(alg.det(S_est)) + np.trace(np.dot(S_est, empCov_set[nodeID])) + K)
    e4.append(alg.norm(S_est -  S_previous, 'fro'))
    
    display = False
    if display == True:
        if (nodeID >timeShift -10) and (nodeID < timeShift + 10):
            print 'nodeID = ', nodeID
            print 'S_true = ', S_actual,'\nS_est', S_est
#            print 'S_error = ',S_actual - S_est, '\n its Fro error = ', alg.norm(S_actual - S_est, 'fro')
            print 'D = ',D,'T = ', T,'TandD = ', TandD,'K = ', K,'P = ', P,'R = ', R,'Score = ', 2* P*R/(P+R)
            
    return e1, e2, e3, e4

def subgradientPolyakCFM(f, sgf, x0, lowerBounder, gamma=1, maxIters=100, report=None):
    ''' Use the Polyak step size and CFM direction update rule. '''
    xk = x0
    gk = sgf(x0)
    sk = gk
    lb = lowerBounder(xk)
    fb = f(xk)
    # Using optimal step size for fixed number of iterations
    for k in arange(maxIters) + 1:
        if report:
            report(xk, k)
        sko = sk
        gko = gk
        gk = sgf(xk)
        #betak = 0.25
        '''betak = (- gamma * gk.dot(sko)/(norm(sko)**2)
                 if sko.dot(gk) < 0 else 0 )'''
        betak = (- gamma * gk.dot(gko) / (norm(gko) ** 2)
                 if gko.dot(gk) < 0 else 0 )
        #sk =  gk + betak*sko
        sk = gk + betak * gko
        nfb = f(xk)
        fb = fb if fb < nfb else nfb
        if nfb < fb:
            nlb = lowerBounder(xk)
            lb = lb if lb > nlb else nlb

        alphak = 0.5 * (fb - lb) / (norm(sk)) #**2 ) #* (k ** -0.66)
        xk = xk - alphak * sk

    return xk

def problem8():
   "problem set 2.1, problem 8, page 56"
   import LUdecomp
   A = np.array([[-3,6,-4],[9,-8,24],[-12,24,-26]],dtype=float)
   A_orig = A.copy()
   LU = LUdecomp.LUdecomp(A)
   b = np.array([-3,65,-42],dtype=float)
   b_orig = b.copy()
   x = LUdecomp.LUsolve(LU,b)
   # extract L and U for verification
   U = np.triu(LU)  # 
   L = np.tril(LU)
   L[ np.diag_indices_from(L) ] = 1.0 
   print("""
Problem 8:
A = 
{}
LU decomposition A = LU, LU (in one matrix) = 
{}
Solving Ax=b, with b = {}
Solution x = {}
Verifying solution: 
     residual ||Ax-b||_2 = {}
     ||A - dot(L,U)||_inf = {}
""".format(A_orig,LU,b_orig,x, 
   la.norm(np.dot(A_orig,x)-b_orig,2), 
   la.norm(A_orig - np.dot(L,U),np.inf))
   )