def test_2d_2d(self):
    #Not dot with vector exactly,
    #just to make sure new feature hasn't break anything

    # Test with row > col
    av = expr.arange((132, 100))
    bv = expr.arange((100, 77))
    na = np.arange(13200).reshape(132, 100)
    nb = np.arange(7700).reshape(100, 77)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

    # Test with row < col
    av = expr.arange((67, 100))
    bv = expr.arange((100, 77))
    na = np.arange(6700).reshape(67, 100)
    nb = np.arange(7700).reshape(100, 77)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

    #Dot with numpy obj
    cv = expr.arange((77, 100))
    dv = np.arange(8800).reshape(100, 88)
    nc = np.arange(7700).reshape(77, 100)
    nd = np.arange(8800).reshape(100, 88)

    Assert.all_eq(expr.dot(cv, dv).glom(),
                  np.dot(nc, nd))

  def test_reshape_dot(self):
    npa1 = np.random.random((357, 93))
    npa2 = np.random.random((31, 357))
    result = np.dot(np.reshape(npa1, (1071, 31)), npa2)

    t1 = expr.from_numpy(npa1)
    t2 = expr.from_numpy(npa2)
    t3 = expr.dot(expr.reshape(t1, (1071, 31)), t2)
    Assert.all_eq(result, t3.glom(), 10e-9)

    npa1 = np.random.random((357, 718))
    npa2 = np.random.random((718, ))
    result = np.dot(npa1, np.reshape(npa2, (718, 1)))

    t1 = expr.from_numpy(npa1)
    t2 = expr.from_numpy(npa2)
    t3 = expr.dot(t1, expr.reshape(t2, (718, 1)))
    Assert.all_eq(result, t3.glom(), 10e-9)

    npa1 = np.random.random((718, ))
    npa2 = np.random.random((1, 357))
    result = np.dot(np.reshape(npa1, (718, 1)), npa2)

    t1 = expr.from_numpy(npa1)
    t2 = expr.from_numpy(npa2)
    t3 = expr.dot(expr.reshape(t1, (718, 1)), t2)
    Assert.all_eq(result, t3.glom(), 10e-9)

def fn2():
  a = expr.ones((N, N))
  b = expr.ones((N, N/2))
  g = expr.dot(a, b) + expr.dot(expr.sum(a, axis=1).reshape((1, N)), b)
  t1 = time.time()
  g_opt = g.optimized()
  #g_opt.force()
  t2 = time.time()
  print t2 - t1
  print g_opt

 def update(self):
     """
 gradient_update = 2xTxw - 2xTy + 2* lambda * w
 Correct this if the update function is wrong.
 """
     xT = expr.transpose(self.x)
     g1 = expr.dot(expr.dot(xT, self.x), self.w)
     g2 = expr.dot(xT, self.y)
     g3 = self.ridge_lambda * self.w
     g4 = g1 + g2 + g3
     return expr.reshape(g4, (1, self.N_DIM))

def bfs(ctx, dim):
	util.log_info("start to computing......")

	sGenerate = time.time()
	current = eager(
			expr.shuffle(
				expr.ndarray(
					(dim, 1),
					dtype = np.int64,
					tile_hint = (dim / ctx.num_workers, 1)),
				make_current,
			))
	
	linkMatrix = eager(
				expr.shuffle(
					expr.ndarray(
					 (dim, dim),
					 dtype = np.int64,
					 tile_hint = (dim, dim / ctx.num_workers)),
				make_matrix,
				))
	eGenerate = time.time()

	startCompute = time.time()
	while(True):
		next = expr.dot(linkMatrix, current)
		formerNum = expr.count_nonzero(current)
		laterNum = expr.count_nonzero(next)
		hasNew = expr.equal(formerNum, laterNum).glom()
		current = next
		if (hasNew):
			break
	current.evaluate()
	endCompute = time.time()
	return (eGenerate - sGenerate, endCompute - startCompute) 

  def _test_optimization_ordered(self):
    na = np.random.rand(1000, 1000)
    nb = np.random.rand(1000, 1000)
    a = expr.from_numpy(na)
    b = expr.from_numpy(nb)

    c = a - b
    d = a + c
    f = c[200:900, 200:900]
    g = d[200:900, 200:900]
    h = f - g
    i = f + h
    j = h[100:500, 100:500]
    k = i[100:500, 100:500]
    l = expr.dot(j, k)
    m = j + k
    n = k - l
    o = n - m
    q = o[100:200, 100:200]

    nc = na - nb
    nd = na + nc
    nf = nc[200:900, 200:900]
    ng = nd[200:900, 200:900]
    nh = nf - ng
    ni = nf + nh
    nj = nh[100:500, 100:500]
    nk = ni[100:500, 100:500]
    nl = np.dot(nj, nk)
    nm = nj + nk
    nn = nk - nl
    no = nn - nm
    nq = no[100:200, 100:200]

    Assert.all_eq(nq, q.optimized().glom(), tolerance = 1e-10)

 def update(self):
   '''
   gradient_update = (h(w) - y) * x
   h(w) = x * w
   '''
   yp = expr.dot(self.x, self.w)
   return self.x * (yp - self.y)

def jacobi_method(A, b, _iter=100):
    """
  Iterative algorithm for approximating the solutions of a diagonally dominant system of linear equations. 

  Parameters
  ----------
  A : ndarray or Expr - 2d
      Input matrix
  b : ndarray or Expr - vector
      RHS vector
  _iter : int
      Times of iteration needed, default to be 100

 Returns
  -------
  result : Expr - vector
      Approximated solution.
  """
    util.Assert.eq(A.shape[0], b.shape[0])

    x = expr.zeros((A.shape[0],))

    D = expr.diag(A)
    R = A - expr.diagflat(D)

    for i in xrange(_iter):
        x = (b - expr.dot(R, x)) / D

    return x

def fn2():
  a = expr.ones((N, N))
  b = expr.ones((N, N))
  x = expr.dot(a, b)
  g = a + b + x

  return g

def cgit(A, x):
  '''
  CGIT Conjugate Gradient iteration
  z = cgit(A, x) generates approximate solution to A*z = x.
  
  Args:
  A(Expr): matrix to be processed.
  x(Expr): the input vector.
  '''
  z = expr.zeros(x.shape, tile_hint=(A.tile_shape()[1], 1))
  r = x
  rho = expr.sum(r * r).glom()
  #util.log_warn('rho:%s', rho)
  p = r
  
  for i in xrange(15):
    q = expr.dot(A, p, tile_hint=(A.tile_shape()[1], 1))
    alpha = rho / expr.sum(p * q).glom()
    #util.log_warn('alpha:%s', alpha)
    z = z + p * alpha
    rho0 = rho
    r = r - q * alpha
    rho = expr.sum(r * r).glom()
    beta = rho / rho0
    #util.log_warn('beta:%s', beta)
    p = r + p * beta
  
  return z

def connectedConponents(ctx, dim, numIters):
	linkMatrix = eager(
					expr.shuffle(
						expr.ndarray(
							(dim, dim),
							dtype = np.int64,
							tile_hint = (dim / ctx.num_workers, dim)),
						make_matrix,
					))

	power = eager(
					expr.shuffle(
						expr.ndarray(
							(dim, dim),
							dtype = np.int64,
							tile_hint = (dim / ctx.num_workers, dim)),
						make_matrix,
					))

	eye = expr.eye(dim, tile_hint = (dim / ctx.num_workers,dim))
	startCompute = time.time()
	result = expr.logical_or(eye, linkMatrix).optimized().glom()
	for i in range(numIters):
		power = expr.dot(power, linkMatrix).optimized().glom()
		result = expr.logical_or(result, power)
	result.optimized().glom()
	final = expr.logical_and(result, expr.transpose(result.optimized())).optimized().evaluate()
	endCompute = time.time()
	return endCompute - startCompute

  def test_numpy_vec_2d(self):
    av = expr.arange(stop = 100)
    bv = np.arange(7700).reshape(100, 77)
    na = np.arange(100)
    nb = np.arange(7700).reshape(100, 77)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

def gen_dot(a, b):
  if not hasattr(a, 'shape') or not hasattr(b, 'shape') or len(a.shape) * len(b.shape) == 0: return [a * b]

  if a.shape[0] == b.shape[0]:
    if len(a.shape) > 1: return [expr.dot(expr.transpose(a), b)]
    elif len(b.shape) == 1: return [expr.dot(a, b)]

  if len(a.shape) > 1 and a.shape[1] == b.shape[0]:
      return [expr.dot(a, b)]

  if len(b.shape) > 1 and a.shape[0] == b.shape[1]:
      return [expr.dot(b, a)]

  if len(a.shape) > 1 and len(b.shape) > 1 and a.shape[1] == b.shape[1]:
      return [expr.dot(a, expr.transpose(b))]

  return [a, b]

 def update(self):
   '''
   gradient_update = (h(w) - y) * x
   h(w) = 1 / (1 + e^(-(x*w)))
   '''
   g = expr.exp(expr.dot(self.x, self.w))
   yp = g / (g + 1)
   return self.x * (yp - self.y)

  def test_numpy_2d_vec(self):
    av = expr.arange((77, 100))
    bv = np.arange(100)
    na = np.arange(7700).reshape(77, 100)
    nb = np.arange(100)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

 def _step():
   yp = expr.dot(x, w)
   Assert.all_eq(yp.shape, y.shape)
   
   diff = x * (yp - y)
   grad = expr.sum(diff, axis=0).glom().reshape((N_DIM, 1))
   wprime = w - grad * 1e-6
   expr.force(wprime)

  def test_numpy_vec_vec(self):
    av = expr.arange(stop=100)
    bv = np.arange(100)
    na = np.arange(100)
    nb = np.arange(100)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

  def train_smo_1998(self, data, labels):
    '''
    Train an SVM model using the SMO (1998) algorithm.
   
    Args:
      data(Expr): points to be trained
      labels(Expr): the correct labels of the training data
    '''
    
    N = data.shape[0] # Number of instances
    D = data.shape[1]  # Number of features

    self.b = 0.0
    self.alpha = expr.zeros((N,1), dtype=np.float64, tile_hint=[N/self.ctx.num_workers, 1]).force()
    
    # linear kernel
    kernel_results = expr.dot(data, expr.transpose(data), tile_hint=[N/self.ctx.num_workers, N])   
    
    labels = expr.force(labels)
    self.E = expr.zeros((N,1), dtype=np.float64, tile_hint=[N/self.ctx.num_workers, 1]).force()
    for i in xrange(N):
      self.E[i, 0] = self.b + expr.reduce(self.alpha, axis=None, dtype_fn=lambda input: input.dtype,
                                          local_reduce_fn=margin_mapper,
                                          accumulate_fn=np.add, 
                                          fn_kw=dict(label=labels, data=kernel_results[:,i].force())).glom() - labels[i, 0]
    
    util.log_info("Starting SMO")
    it = 0
    num_changed = 0
    examine_all = True
    while (num_changed > 0 or examine_all) and (it < self.maxiter):
      util.log_info("Iteration:%d", it)

      num_changed = 0
      
      if examine_all:
        for i in xrange(N): 
          num_changed += self.examine_example(i, N, labels, kernel_results)
      else:
        for i in xrange(N):
          if self.alpha[i, 0] > 0 and self.alpha[i, 0] < self.C:
            num_changed += self.examine_example(i, N, labels, kernel_results)

      it += 1

      if examine_all: examine_all = False
      elif num_changed == 0: examine_all = True
    
    self.w = expr.zeros((D, 1), dtype=np.float64).force()
    for i in xrange(D): 
      self.w[i,0] = expr.reduce(self.alpha, axis=None, dtype_fn=lambda input: input.dtype,
                              local_reduce_fn=margin_mapper,
                              accumulate_fn=np.add, 
                              fn_kw=dict(label=labels, data=expr.force(data[:,i]))).glom()
    self.usew_ = True
    print 'iteration finish:', it
    print 'b:', self.b
    print 'w:', self.w.glom()

def svd(A, k=None):
  """
  Stochastic SVD.

  Parameters
  ----------
  A : spartan matrix
      Array to compute the SVD on, of shape (M, N)
  k : int, optional
      Number of singular values and vectors to compute.

  The operations include matrix multiplication and QR decomposition.
  We parallelize both of them.

  Returns
  --------
  U : Spartan array of shape (M, k)
  S : numpy array of shape (k,)
  V : numpy array of shape (k, k)
  """
  if k is None:
    k = A.shape[1]

  ctx = blob_ctx.get()
  Omega = expr.randn(A.shape[1], k, tile_hint=(A.shape[1]/ctx.num_workers, k))

  r = A.shape[0] / ctx.num_workers
  Y = expr.dot(A, Omega, tile_hint=(r, k)).force()
  
  Q, R = qr(Y)
  
  B = expr.dot(expr.transpose(Q), A)
  BTB = expr.dot(B, expr.transpose(B)).glom()

  S, U_ = np.linalg.eig(BTB)
  S = np.sqrt(S)

  # Sort by eigen values from large to small
  si = np.argsort(S)[::-1]
  S = S[si]
  U_ = U_[:, si]

  U = expr.dot(Q, U_).force()
  V = np.dot(np.dot(expr.transpose(B).glom(), U_), np.diag(np.ones(S.shape[0]) / S))
  return U, S, V.T 

  def test_2d_vec(self):
    # Test with row > col
    av = expr.arange((100, 77))
    bv = expr.arange(stop=77)
    na = np.arange(7700).reshape(100, 77)
    nb = np.arange(77)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))

    # Test with col > row
    av = expr.arange((77, 100))
    bv = expr.arange(stop=100)
    na = np.arange(7700).reshape(77, 100)
    nb = np.arange(100)

    Assert.all_eq(expr.dot(av, bv).glom(),
                  np.dot(na, nb))