def fmm_single_wall_stokeslet(r_vectors, force, eta, a, *args, **kwargs):
  '''
  WARNING: pseudo-PBC are not implemented for this function.

  Compute the Stokeslet interaction plus self mobility
  II/(6*pi*eta*a) in the presence of a wall at z=0.
  It uses the fmm implemented in the library stfmm3d.
  Must compile mobility_fmm.f90 before this will work
  (see Makefile).

  For blobs overlaping the wall we use
  Compute M = B^T * M_tilde(z_effective) * B.
  '''
  # Get effective height
  r_vectors_effective = shift_heights(r_vectors, a)
  # Compute damping matrix B
  B, overlap = damping_matrix_B(r_vectors, a, *args, **kwargs)
  # Compute B * force
  if overlap is True:
    force = B.dot(force)
  # Compute M_tilde * B * vector
  num_particles = r_vectors.size // 3
  ier = 0
  iprec = 5
  r_vectors_fortran = np.copy(r_vectors_effective.T, order='F')
  force_fortran = np.copy(np.reshape(force, (num_particles, 3)).T, order='F')
  u_fortran = np.empty_like(r_vectors_fortran, order='F')
  fmm.fmm_stokeslet_half(r_vectors_fortran, force_fortran, u_fortran, ier, iprec, a, eta, num_particles)
  # Compute B.T * M * B * force
  if overlap is True:
    return B.dot(np.reshape(u_fortran.T, u_fortran.size))
  else:
    return np.reshape(u_fortran.T, u_fortran.size)

 def make_net(self, input_images, input_measurements, input_actions, input_objectives, reuse=False):
     if reuse:
         tf.get_variable_scope().reuse_variables()
     
     self.fc_val_params = np.copy(self.fc_joint_params)
     self.fc_val_params['out_dims'][-1] = self.target_dim
     self.fc_adv_params = np.copy(self.fc_joint_params)
     self.fc_adv_params['out_dims'][-1] = len(self.net_discrete_actions) * self.target_dim
     p_img_conv = my_ops.conv_encoder(input_images, self.conv_params, 'p_img_conv', msra_coeff=0.9)
     p_img_fc = my_ops.fc_net(my_ops.flatten(p_img_conv), self.fc_img_params, 'p_img_fc', msra_coeff=0.9)
     p_meas_fc = my_ops.fc_net(input_measurements, self.fc_meas_params, 'p_meas_fc', msra_coeff=0.9)
     if isinstance(self.fc_obj_params, np.ndarray):
         p_obj_fc = my_ops.fc_net(input_objectives, self.fc_obj_params, 'p_obj_fc', msra_coeff=0.9)
         p_concat_fc = tf.concat([p_img_fc,p_meas_fc,p_obj_fc], 1)
     else:
         p_concat_fc = tf.concat([p_img_fc,p_meas_fc], 1)
         if self.random_objective_coeffs:
             raise Exception('Need fc_obj_params with randomized objectives')
         
     p_val_fc = my_ops.fc_net(p_concat_fc, self.fc_val_params, 'p_val_fc', last_linear=True, msra_coeff=0.9)
     p_adv_fc = my_ops.fc_net(p_concat_fc, self.fc_adv_params, 'p_adv_fc', last_linear=True, msra_coeff=0.9)
     
     adv_reshape = tf.reshape(p_adv_fc, [-1, len(self.net_discrete_actions), self.target_dim])
     
     pred_all_nomean = adv_reshape - tf.reduce_mean(adv_reshape, reduction_indices=1, keep_dims=True)
     pred_all = pred_all_nomean + tf.reshape(p_val_fc, [-1, 1, self.target_dim])
     pred_relevant = tf.boolean_mask(pred_all, tf.cast(input_actions, tf.bool))
     
     return pred_all, pred_relevant

def replica_exchange(index_spin, index_replica, beta, d_beta,X = [[]] ):
    x1, x2 = X[index_replica, :], X[index_replica + 1, :]
    # r = P(x_k | beta_k+1)P(x_k+1 | beta_k) / P(x_k | beta_k)P(x_k+1 | beta_k+1)
    r = beta_power_of_prob(index_spin, beta + d_beta, x1, theta) * beta_power_of_prob(index_spin, beta, x2, theta) / beta_power_of_prob(index_spin, beta, x1, theta) * beta_power_of_prob(index_spin, beta +d_beta, x2, theta)
    if(np.random.uniform(size=1) < r):
        X[index_replica, :], X[index_replica, :] = np.copy(x2), np.copy(x1)
    return X

 def clear(self):
     self.new_x = copy(self.x)
     self.new_y = copy(self.y)
     args = ()
     for i in range(len(self.x)):
         args = args + (self.new_x[i], self.new_y[i])
     self.plot(*args)

def copyCase(case):
    ppc = {"version": 2}
    ppc["baseMVA"] = 100.0
    ppc["bus"] = copy(case["bus"])
    ppc["gen"] = copy(case["gen"])
    ppc["branch"] = copy(case["branch"])
    return ppc;

    def make_corners(self, f):
        """
        The standard mom grid includes t-cell corners be specifying the u, v
        grid. Here we extract that and put it into the format expected by
        the regridder and OASIS.
        """

        x = np.copy(f.variables['x'])
        y = np.copy(f.variables['y'])

        self.clon = np.zeros((4, x.shape[0] / 2, x.shape[1] / 2))
        self.clon[:] = np.NAN
        self.clat = np.zeros((4, x.shape[0] / 2, x.shape[1] / 2))
        self.clat[:] = np.NAN

        # Corner lats. 0 is bottom left and then counter-clockwise. 
        # This is the OASIS convention. 
        self.clat[0,:,:] = y[0:-1:2,0:-1:2]
        self.clat[1,:,:] = y[0:-1:2,2::2]
        self.clat[2,:,:] = y[2::2,2::2]
        self.clat[3,:,:] = y[2::2,0:-1:2]

        # Corner lons.
        self.clon[0,:,:] = x[0:-1:2,0:-1:2]
        self.clon[1,:,:] = x[0:-1:2,2::2]
        self.clon[2,:,:] = x[2::2,2::2]
        self.clon[3,:,:] = x[2::2,0:-1:2]

        # Select points from double density grid. Southern most U points are
        # excluded. also the last (Eastern) U points, they are duplicates of
        # the first.
        self.x_t = x[1::2,1::2]
        self.y_t = y[1::2,1::2]
        self.x_u = x[2::2,0:-1:2]
        self.y_u = y[2::2,0:-1:2]

def lossFun(inputs, targets, hprev):
  """
  inputs,targets are both list of integers.
  hprev is Hx1 array of initial hidden state
  returns the loss, gradients on model parameters, and last hidden state
  """
  xs, hs, ys, ps = {}, {}, {}, {}
  hs[-1] = np.copy(hprev)
  loss = 0
  # forward pass
  for t in xrange(len(inputs)):
    xs[t] = np.zeros((vocab_size,1)) # encode in 1-of-k representation
    xs[t][inputs[t]] = 1
    hs[t] = np.tanh(np.dot(Wxh, xs[t]) + np.dot(Whh, hs[t-1]) + bh) # hidden state
    ys[t] = np.dot(Why, hs[t]) + by # unnormalized log probabilities for next chars
    ps[t] = np.exp(ys[t]) / np.sum(np.exp(ys[t])) # probabilities for next chars
    loss += -np.log(ps[t][targets[t],0]) # softmax (cross-entropy loss)
  # backward pass: compute gradients going backwards
  dWxh, dWhh, dWhy = np.zeros_like(Wxh), np.zeros_like(Whh), np.zeros_like(Why)
  dbh, dby = np.zeros_like(bh), np.zeros_like(by)
  dhnext = np.zeros_like(hs[0])
  for t in reversed(xrange(len(inputs))):
    dy = np.copy(ps[t])
    dy[targets[t]] -= 1 # backprop into y
    dWhy += np.dot(dy, hs[t].T)
    dby += dy
    dh = np.dot(Why.T, dy) + dhnext # backprop into h
    dhraw = (1 - hs[t] * hs[t]) * dh # backprop through tanh nonlinearity
    dbh += dhraw
    dWxh += np.dot(dhraw, xs[t].T)
    dWhh += np.dot(dhraw, hs[t-1].T)
    dhnext = np.dot(Whh.T, dhraw)
  for dparam in [dWxh, dWhh, dWhy, dbh, dby]:
    np.clip(dparam, -5, 5, out=dparam) # clip to mitigate exploding gradients
  return loss, dWxh, dWhh, dWhy, dbh, dby, hs[len(inputs)-1]

def classify(image, hog, rho, max_detected=8):
    image_boxes = np.copy(image)
    found = hog.detect(image_boxes, winStride=(1, 1))

    if len(found[0]) == 0:
        return "female", image_boxes, 0

    scores = np.zeros(found[1].shape[0])
    for index, score in enumerate(found[1]):
        scores[index] = found[1][index][0]
    order = np.argsort(scores)

    image_boxes = np.copy(image)
    index = 0
    while index < max_detected and found[1][order[index]] - rho < 0:
        current = found[0][order[index], :]
        x, y = current
        h = hog.compute(image[y : (y + win_height), x : (x + win_width), :])
        colour = (0, 255, 0)
        cv2.rectangle(image_boxes, (x, y), (x + win_width, y + win_height), colour, 1)
        index += 1
        # print 'Number of detected objects = %d' % index

    return (
        "male" if index > 0 else "female",
        image_boxes,
        index,
        found[0][order[(index - 1) : index], :],
        found[1][order[(index - 1) : index]],
    )

def get_centroids(points, k):
    '''KMeans++的初始化聚类中心的方法
    input:  points(mat):样本
            k(int):聚类中心的个数
    output: cluster_centers(mat):初始化后的聚类中心
    '''
    m, n = np.shape(points)
    cluster_centers = np.mat(np.zeros((k , n)))
    # 1、随机选择一个样本点为第一个聚类中心
    index = np.random.randint(0, m)
    cluster_centers[0, ] = np.copy(points[index, ])
    # 2、初始化一个距离的序列
    d = [0.0 for _ in xrange(m)]
 
    for i in xrange(1, k):
        sum_all = 0
        for j in xrange(m):
            # 3、对每一个样本找到最近的聚类中心点
            d[j] = nearest(points[j, ], cluster_centers[0:i, ])
            # 4、将所有的最短距离相加
            sum_all += d[j]
        # 5、取得sum_all之间的随机值
        sum_all *= random()
        # 6、获得距离最远的样本点作为聚类中心点
        for j, di in enumerate(d):
            sum_all -= di
            if sum_all > 0:
                continue
            cluster_centers[i] = np.copy(points[j, ])
            break
    return cluster_centers

def Solver2(A,m,n,disc,Sol,score):
    "Burst the component with maximum bubbles"
    score=0
    Acopy=np.copy(A)#copy of A.So that changes to A aren't reflected outside.Copy by value
    #I=0
    while(1):
        #I=I+1
        #print str(I)+"th iteration"
        components=[]
        disc[:]=0
        for i in xrange(m):
            for j in xrange(n):
                if Acopy[i][j]!=-1:# an unburst bubble
                    l=bfs(Acopy,i,j,disc,m,n)
                    components.append(l)
        #print "No of components=",len(components)            
        if len(components)>0:            
            maxim=len(components[0])
        else:#game over all bubble gone
            #print "Game over all finished."
            Sol.append(score)
            break    
        b=components[0]#initializing component which will be burst
        for c in components:
            if len(c)>maxim:
                maxim=len(c)
                b=c
        if maxim==1:#this means game over
            #print "Game over"
            Sol.append(score)
            break
        burst(Acopy,b,m,n)
        Mcopy=np.copy(Acopy)
        Sol.append([Mcopy,len(b)**2])
        score=score+len(b)**2            

def Solver3(A,m,n,disc,Sol,score):
    "Burst a component at random"
    score=0
    Acopy=np.copy(A)
    while(1):
        #I=I+1
        #print str(I)+"th iteration"
        components=[]
        disc[:]=0
        for i in xrange(m):
            for j in xrange(n):
                if Acopy[i][j]!=-1:# an unburst bubble
                    l=bfs(Acopy,i,j,disc,m,n)
                    components.append(l)
        #print "No of components=",len(components)            
        if len(components)==0:            
            Sol.append(score)
            break
        breakable=[]#contains those components which can be burst    
        for c in components:
            if len(c)>1:
               breakable.append(c)
        if len(breakable)==0:
            Sol.append(score)
            break               
        b=breakable[random.randrange(0,len(breakable))]#random component which will be burst
        burst(Acopy,b,m,n)
        Mcopy=np.copy(Acopy)
        Sol.append([Mcopy,len(b)**2])
        score=score+len(b)**2       

def local_search( start_node, goal_cost=0 ):
    node = start_node
    
    while node.cost > goal_cost:
      
      while True:
         """
         Pick a random row and check if it is causing conflicts
         """
         i = random.randint( 0, start_node.gene.shape[0]-1 ) # select row to manipulate
         tmp = np.copy( node.gene )
         tmp[i,:]=0
         if objective_function(tmp)<node.cost:
            # Yep, this caused some conflicts
            break
      
      neighbors = [ (float("inf"),None) ]

      for modified in generate_permutations(node.gene[i]):
        tmp = np.copy( node.gene )
        tmp[i] = modified
        node = Board(tmp)
        
        if node.cost == neighbors[0][0]:
           neighbors.append((node.cost,node))
        elif node.cost<neighbors[0][0]:
           neighbors = [ (node.cost,node)]
           
      
      node = random.choice(neighbors)[1]
    #end while
    
    return node

def hun(costMatrix):

    # Check first, if costmatrix is not empty
    if costMatrix.shape==(0,0):
        return []

    # Create squared temporary matrix
    tmpMatrix = numpy.copy(costMatrix)
    tmpMatrix = makeSquareWithNegValues(tmpMatrix)
    sqCostMatrix = numpy.copy(tmpMatrix)
    sqCostMatrix[tmpMatrix==-1]=10e10

    # Solve ASP on the temporary matrix
    m=Munkres()
    i=m.compute(sqCostMatrix)


    # Create resultin matrix that contains ones at matching
    # objects and remove all excluded matches
    binMatrix = numpy.zeros( tmpMatrix.shape,dtype=bool )
    for x,y in i:
        if tmpMatrix[x,y]==-1:
            continue
        binMatrix[x,y]=True

    return binMatrix

    def getCurrentSpectrum(self):
        # self.c_rfi.execute("SELECT spectrum, timestamp from spectra_%i where timestamp = (SELECT max(timestamp) from spectra_%i)"%(self.which_db, self.which_db))
        # result = self.c_rfi.fetchone()
        # self.c_rfi.execute("FLUSH QUERY CACHE")

        data = numpy.copy(self.curr)

        result = [cnf.remove_internal_RFI(data,self.mode[0]), self.time[0], self.mode[0]]

        # print "current timestamp = %i"%self.time[0]

        while self.last_timestamp == result[1]: #Current timestamp equals last timestamp
            data = numpy.copy(self.curr)
            result = [cnf.remove_internal_RFI(data,self.mode[0]), self.time[0], self.mode[0]]
            time.sleep(0.1)
            # res = self.c_rfi.fetchall()
            # if res[0][0] != self.which_db: #Check if using the correct DB
            #     self.which_db = res[0][0]
            # while self.last_timestamp == result[1]:
            #     self.c_rfi.execute("SELECT spectrum, timestamp from spectra_%i where timestamp = (SELECT max(timestamp) from spectra_%i)"%(self.which_db, self.which_db))
            #     result = self.c_rfi.fetchone()
            #     self.c_rfi.execute("FLUSH QUERY CACHE")
            #     time.sleep(0.1)
        self.last_timestamp = result[1]
        return (result[0], result[1], result[2] - 1)

    def cluster(self, data, n_clusters):

        n, d = shape(data)
        locations = zeros((self.n_particles, n_clusters, d))

        for i in range(self.n_particles):
            for j in range(n_clusters):
                locations[i, j, :] = copy(data[randint(n), :])  # Initialize cluster centers to random datapoints

        bestlocations = copy(locations)
        velocities = zeros((self.n_particles, n_clusters, d))

        bestscores = [score(data, centroids=locations[i, :, :], norm=self.norm) for i in range(self.n_particles)]
        sbestlocation = copy(locations[argmin(bestscores), :, :])
        sbestscore = min(bestscores)

        for i in range(self.n_iterations):
            if i % self.printfreq == 0:
                print "Particle swarm iteration", i, "best score:", sbestscore
            for j in range(self.n_particles):
                r = rand(n_clusters, d)
                s = rand(n_clusters, d)
                velocities[j, :, :] = (self.w * velocities[j, :, :]) + \
                                      (self.c1 * r * (bestlocations[j, :, :] - locations[j, :, :])) + \
                                      (self.c2 * s * (sbestlocation - locations[j, :, :]))
                locations[j, :, :] = locations[j, :, :] + velocities[j, :, :]
                currentscore = score(data, centroids=locations[j, :, :], norm=self.norm)
                if currentscore < bestscores[j]:
                    bestscores[j] = currentscore
                    bestlocations[j, :, :] = locations[j, :, :]
                    if currentscore < sbestscore:
                        sbestscore = currentscore
                        sbestlocation = copy(locations[j, :, :])

        return getlabels(data, centroids=sbestlocation, norm=self.norm)

 def get_unidirectional_S(self):
     S_plus = np.copy(self.S)
     S_minus = np.copy(self.S)
     S_plus[self.S < 0] = 0
     S_minus[self.S > 0] = 0
     return S_minus, S_plus
     

 def __array__(self, dtype=None):
     if self.size:
         arrayfire.backend.get().af_get_data_ptr(ctypes.c_void_p(self.h_array.ctypes.data), self.d_array.arr)
     if dtype is None:
         return numpy.copy(self.h_array)
     else:
         return numpy.copy(self.h_array).astype(dtype)

def cap(guess_vector):
    """
    This takes the Euler equations, and sets them equal to zero for an f-solve
    Remember that Keq was found by taking the derivative of the sum of the 
        utility functions, with respect to k in each time period, and that 
        leq was the same, but because l only shows up in 1 period, it has a
        much smaller term.

    ### Paramaters ###
    guess_vector: The first half is the intial guess for the kapital, and
        the second half is the intial guess for the labor
    """
    #equations for keq
    ks = np.zeros(periods)
    ks[1:] = guess_vector[:periods-1]
    ls  = guess_vector[periods-1:]
    kk  = ks[:-1]
    kk1 = ks[1:]
    kk2 = np.zeros(periods-1)
    kk2[:-1] = ks[2:]
    lk  = ls[:-1]
    lk1 = ls[1:]
    #equation for leq
    ll = np.copy(ls)
    kl = np.copy(ks)
    kl1 = np.zeros(periods)
    kl1[:-1] = kl[1:]
    w = wage(ks, ls)
    r = rate(ks, ls)
    keq = ((lk*w+(1.+r-delta)*kk - kk1)**-gamma - (beta*(1+r-delta)*(lk1*w+(1+r-delta)*kk1-kk2)**-gamma))
    leq = ((w*(ll*w + (1+r-delta)*kl-kl1)**-gamma)-(1-ll)**-sigma)
    error = np.append(keq, leq)

    return np.append(keq, leq)

def build_models2(X):
	mu_init = [-13, -4, 50]
	sigmasq_init = [5, 14, 30] 
	wt_init = [0.2, 0.4, 0.4]
	its = 20

	L = []
	mu = np.copy(mu_init)
	sigmasq = np.copy(sigmasq_init)
	wt = np.copy(wt_init)

	#firt iteration
	result = gmmest(X, mu_init, sigmasq_init, wt_init, its)
	mu = np.array(result[0][:])
	sigmasq = np.array(result[1][:])
	wt = np.array(result[2][:])
	L.append(result[3])

	#rest of iterations
	for i in range(its-1):
		result = gmmest(X, mu, sigmasq, wt, 1)
		mu = np.array(result[0][:])
		sigmasq = np.array(result[1][:])
		wt = np.array(result[2][:])
		L.append(result[3])

	
	return result, L

def _compute_positions(spin, seed_pos, seed_spin):
    """Given an array of SPINS and two SEED_POSITION and SEED_SPIN values, compute the positions along the prime hex 

    """

    logger.info("compute_positions: starting aux calculations")
    delta = np.zeros_like(spin)
    delta[0] = spin[0] - seed_spin  # first delta is cur_spin - prev_spin from seed_spin
    delta[1:] = spin[1:] - spin[:-1]  # delta is cur_spin - prev_spin
    logger.info("compute_positions: delta={}".format(delta))

    increments = np.copy(spin)  # copy the spin array,
    increments[delta != 0] = 0  # set any non-zero delta to zero in the increment array
    logger.info("compute_positions: increments={}".format(increments))

    logger.info("compute_positions:\tdone with aux calculations")

    logger.info("compute_positions: starting primary calculation")

    # start at seed, cumulative add
    positions = np.copy(increments)
    # increments[0] = seed_pos
    outpositions = (seed_pos + np.cumsum(increments)) % 6
    logger.info("compute_positions: outpositions={}".format(outpositions))
    logger.info("compute_positions:\tdone with primary calculation")
    return outpositions