def setdata(self, X, V):
        A = self.bialtprodeye(2*self.F.J_coords)
        """Note: p, q <= min(n,m)"""

        self.data.Brand = 2*(random((A.shape[0],self.data.p))-0.5)
        self.data.Crand = 2*(random((A.shape[1],self.data.q))-0.5)
        self.data.B = zeros((A.shape[0],self.data.p), float)
        self.data.C = zeros((A.shape[1],self.data.q), float)
        self.data.D = zeros((self.data.q,self.data.p), float)

        U, S, Vh = linalg.svd(A)
        self.data.b = U[:,-1:]
        self.data.c = transpose(Vh)[:,-1:]

        if self.update:
            self.data.B[:,1] = self.data.b
            self.data.C[:,1] = self.data.c

            U2, S2, Vh2 = linalg.svd(c_[r_[A, transpose(self.data.C[:,1])], r_[self.data.B[:,1], [[0]]]])
            self.data.B[:,2] = U2[0:A.shape[0],-1:]
            self.data.C[:,2] = transpose(Vh2)[0:A.shape[1],-1:]
            self.data.D[0,1] = U2[A.shape[0],-1]
            self.data.D[1,0] = transpose(Vh2)[A.shape[1],-1]
        else:
                # self.data.B = eye(self.data.Brand.shape)
                # self.data.C = eye(self.data.Crand.shape)
                # USE OF RANDOM
            self.data.B = self.data.Brand
            self.data.C = self.data.Crand

 def _testSpeednD(self,nr_points, dim):
     points = 100*numpy.random(dim, nr_points)
     vals = 100*numpy.random(nr_points)
     
     #build the model
     lut_ss = LutStrategy()
     lut_model = LutFactory().build(points, vals, lut_ss)       
     
     # test the model    
     target_points = points + 0.1
     cnt=0
     
     starttime=time.time()
     
     while cnt < 2:
         yhat = lut_model.simulate(target_points)           
         cnt=cnt+1
     
     elapsed=time.time()-starttime
     
     nb_lookups=nr_points * cnt
     
     lookups_per_sec=nb_lookups / elapsed
     
     print "%d simulations (%d-D) of %d points took %f seconds (%d lookups/sec)" % ( cnt, dim , nr_points, elapsed, lookups_per_sec)

    def __init__(self, center, radius, theshold, maxDepth):

        self.isLeaf = False
        self.center = center
        self.radius = radius
        self.threshold = theshold
        self.maxDepth = maxDepth
        #        self.color = (1., 0., 0.)
        #        self.color = (randint(0,255), randint(0, 255), randint(0, 255))
        self.color = (random(), random(), random())
        # print "self.color = ", self.color

        self.children = [None, None, None, None, None, None, None, None]
        #        self.boundsOffsetTable = array([[-0.5, -0.5, -0.5],
        #                                        [+0.5, -0.5, -0.5],
        #                                        [-0.5, +0.5, -0.5],
        #                                        [+0.5, +0.5, -0.5],
        #                                        [-0.5, -0.5, +0.5],
        #                                        [+0.5, -0.5, +0.5],
        #                                        [-0.5, +0.5, +0.5],
        #                                        [+0.5, +0.5, +0.5]])
        self.boundsOffsetTable = array(
            [
                [-0.5, -0.5, -0.5],
                [+0.5, -0.5, -0.5],
                [-0.5, -0.5, +0.5],
                [+0.5, -0.5, +0.5],
                [-0.5, +0.5, -0.5],
                [+0.5, +0.5, -0.5],
                [-0.5, +0.5, +0.5],
                [+0.5, +0.5, +0.5],
            ]
        )

def update_state():
	newly_infected.a = False
	global has_infection

	# visit the nodes in random order
	vs = list(g.vertices())
	shuffle(vs)

	for v in vs:
		if random() < x:
			p[v] = 0.0
			newly_infected[v] = True
		elif has_infection[v] == True:
			ns = list(v.out_neighbours())
			if len(ns) > 0:
				for w in ns:
					if random() < p[w]: # chance de ser infectado				
						newly_infected[w] = True
						p[w] = 0.0
		if (p[v] < 1.0):
			p[v] += recovery_rate
		state[v] = [p[v], p[v], p[v], 1.0]

	has_infection = cp.deepcopy(newly_infected)	

    #with open("plot.txt", "a") as myfile:
    #	myfile.write(str(counter["S"])+","+str(counter["I"])+","+str(counter["R"])+"\n")
    #plot_values["S"].append(counter["S"])
    #plot_values["I"].append(counter["I"])
    #plot_values["R"].append(counter["R"])

    # Filter out the recovered vertices
    #g.set_vertex_filter(removed, inverted=True)

    # The following will force the re-drawing of the graph, and issue a
    # re-drawing of the GTK window.
	win.graph.regenerate_surface()
	win.graph.queue_draw()


    #ax1.plot(range(len(plot_values["S"])), plot_values["S"],color='b')
    #ax1.plot(range(len(plot_values["I"])), plot_values["I"],color='r')
    #ax1.plot(range(len(plot_values["R"])), plot_values["R"],color='y')
    #fig1.canvas.draw()

    # if doing an offscreen animation, dump frame to disk
	if offscreen:
		global count
		pixbuf = win.get_pixbuf()
		pixbuf.savev(r'./frames/sirs%06d.png' % count, 'png', [], [])
		if count > max_count:
			sys.exit(0)
		count += 1

    # We need to return True so that the main loop will call this function more
    # than once.
	return True

def mutate(p):
    if random() < 0.5:
        return random(), random()
    else:
        dx = normal(0, 0.0001)
        dy = normal(0, 0.0001)
        x = p[0] + dx
        y = p[1] + dy
        x = x - 1 if x > 1 else x + 1 if x < 0 else x
        y = y - 1 if y > 1 else y + 1 if y < 0 else y
        return x, y

def update_state():
    newly_infected.a = False
    removed.a = False

    # visit the nodes in random order
    vs = list(g.vertices())
    shuffle(vs)
    for v in vs:
        if state[v] == I:
            if random() < r:
                state[v] = R
        elif state[v] == S:
            if random() < x:
                state[v] = I
            else:
                ns = list(v.out_neighbours())
                if len(ns) > 0:
                    w = ns[randint(0, len(ns))]  # choose a random neighbour
                    if state[w] == I:
                        state[v] = I
                        newly_infected[v] = True
        elif random() < s:
            state[v] = S
        if state[v] == R:
            removed[v] = True

        if state[v] == S:
            if I in [state[w] for w in v.out_neighbours()]:
                vertex_sfcs[v] = Simg_fear
            else:
                vertex_sfcs[v] = Simg
        else:
            vertex_sfcs[v] = Iimg

    # Filter out the recovered vertices
    g.set_vertex_filter(removed, inverted=True)

    # The following will force the re-drawing of the graph, and issue a
    # re-drawing of the GTK window.
    win.graph.regenerate_surface(lazy=False)
    win.graph.queue_draw()

    # if doing an offscreen animation, dump frame to disk
    if offscreen:
        global count
        pixbuf = win.get_pixbuf()
        pixbuf.savev(r'./frames/zombies%06d.png' % count, 'png', [], [])
        if count > max_count:
            sys.exit(0)
        count += 1

    # We need to return True so that the main loop will call this function more
    # than once.
    return True

 def RandSelectFocalPoint(self):
     objId = random.randint(0, self.GetObjectsCount()-1)
     faceId = random.randint(0, self.GetObjectFacesCount(objId)-1)
     
     pt1, pt2 ,pt3 = self.GetVertices(objId, faceId)
     c1 = numpy.random()
     c2 = numpy.random()
     c3 = numpy.random()
     c = c1 + c2 + c3
     
     pt = (c1/c) * pt1 + (c2/c) * pt2 + (c3/c) * pt3
     return pt

    def setdata(self, A):
        """Note: p, q <= min(n,m)"""
        self.data.Brand = 2*(random((A.shape[0],self.data.p))-0.5)
        self.data.Crand = 2*(random((A.shape[1],self.data.q))-0.5)
        self.data.D = zeros((self.data.q,self.data.p), float)

        if self.update:
            U, S, Vh = linalg.svd(A)
            self.data.B = U[:,-1*self.data.p:]
            self.data.C = transpose(Vh)[:,-1*self.data.q:]
        else:
            self.data.B = self.data.Brand
            self.data.C = self.data.Crand

 def generate_harmonic_oscillators(self, number_ = None):
     if number_ == None:
         self.harmonic_oscillators = \
                 self.x_range*np.random.random(len(self.harmonic_oscillators)) \
                 - self.x_range/2.
     else:
         self.harmonic_oscillators = self.x_range*np.random(number_) - self.x_range/2.

    def step(self):
        # logpability and loglike for stoch's current value:
        logp = sum([stoch.logp for stoch in self.stochs]) + self.indicator.logp
        loglike = self.loglike

        # Sample a candidate value for the value and indicator of the stoch.
        self.propose()

        # logpability and loglike for stoch's proposed value:
        logp_p = sum([stoch.logp for stoch in self.stochs]) + self.indicator.logp

        # Skip the rest if a bad value is proposed
        if logp_p == -Inf:
            for stoch in self.stochs: stoch.revert()
            return

        loglike_p = self.loglike

        # test:
        test_val =  logp_p + loglike_p - logp - loglike
        test_val += self.inv_q(self.indicator)
        test_val += self.q(self.indicator,self._u)

        if self.Jacobian is not None:
            test_val += self.Jacobian(self.indicator,self._u,**self.stoch_dict)

        if log(random()) > test_val:
            for stoch in self.stochs:
                stoch.revert

def hamil_i(numsites, site, J, V, m):
	
	#Params
	#numsites 	- positive integer for the number of sites in spin chain
	#site		- non-negative integer representing the site operated on
	#J,V,h		- arbitrary constants that define system potentials
	#returns (qutip.Qobj) hamiltonian for fermions jumping from site
	#to site oper_i means combined hilbert space operator acting at 
	#site i, with remaining sites untouched by the total operator

	
	#build ith state operators
	sigplus_i = operator_i('sigma plus', numsites, site)
	sigminus_i = operator_i('sigma minus', numsites, site)
	number_i = operator_i('number', numsites, site)
	sigma_zi = operator_i('sigmaz', numsites, site)
	sigplus_ip1 = operator_i('sigma plus', numsites, site+1)
	sigminus_ip1 = operator_i('sigma minus', numsites, site+1)
	number_ip1 = operator_i('number', numsites, site+1)

	#crate parts of the hamiltonian 
	#i.e. H_i = J*(s+_i*s-_1+1 + s-_i*s+_i+1) + V*m_i*m_1+1 + m_i*h_i
	rand_dist = 2*m*(random()-0.5) 
	
	jump = J*(sigplus_i*sigminus_ip1 + sigminus_i*sigplus_ip1)
	interaction_potential = V*(number_i*number_ip1)
	site_potential = rand_dist*sigma_zi #write this into operator_i
	
	#add pieces together
	H_i = jump + interaction_potential + site_potential 
	
	return H_i

def buildArrays( ):
	a = arange(0,n)
	vertex = shuffle(cos(2*pi*a/n), sin(2*pi*a/n))
	vertex.shape = (n, 2)
	color = random(n*3)
	color.shape = (n, 3)
	return vertex,color

    def isample_without_replacement(self, k):
        """ Return a sample of size k, without replacement

        k <= n

        O(n)

        Use a heap to keep track of selection.
        """
        if k > len(self.weights):
            raise ValueError("Sample size should be <= %d" % len(self.weights))
    
        heap = []

        random = self.random.random_sample
        weights = random(len(self.weights)) ** (1.0/self.weights)

        for ix, weight in enumerate(weights):
            if ix < k:
                heapq.heappush(heap, (weight, ix))
            else:
                if heap[0][0] < weight:
                    heapq.heapreplace(heap, (weight, ix))

        # now sort the heap -- this is to make things repeatable
        heap.sort()

        # return permuted indices
        return(self.random.permutation([x[1] for x in heap]))

    def __init__(self, sizes):
        self.sizes = sizes
        self.n_layers = len(sizes) 
        self.learningRate = 2  # Note: typically needs to be lower when using 'sigm' activation function and non-normalized inputs.
        self.momentum = 0.5
        self.scaling_learningRate = 1  # Scaling factor for the learning rate (each epoch)
        self.weightPenaltyL2 = 0 # L2 regularization
        self.nonSparsityPenalty  = 0 #  Non sparsity penalty
        self.sparsityTarget = 0.05 #  Sparsity target
        self.dropoutFraction = 0 # Dropout level
        self.activation_function = 'tanh_opt' # Activation functions of hidden layers: 'sigm' (sigmoid) or 'tanh' (optimal tanh).
        self.output = 'sigm'  # output unit 'sigm' (=logistic), 'softmax' and 'linear'
        self.testing = False    
        self.W = [None for _ in range(1, self.n_layers)]
        self.vW = [None for _ in range(1, self.n_layers)]
        self.p = [None for _ in range(1, self.n_layers)]
        self.n_outputs = self.sizes[-1]

        for i in range(1, self.n_layers):
            # weights and weight momentum
            # +1 in shape for bias
            self.W[i - 1] = (np.random(self.sizes[i], self.sizes[i - 1]+1) - 0.5) * 2 * 4 * sqrt(6 / (self.sizes[i] + self.sizes[i - 1])) 
            self.vW[i - 1] = np.zeros_like(self.W[i - 1])

            # average activations
            self.p[i]= np.zeros(1, self.sizes[i])

	def generatePoint(self, i):
		"""
		Generates a new point through mutation and crossover. 
		"""
		# Select 3 distinct indices different from i between 0 and np-1
		indices=arange(self.Np)
		indices=concatenate((indices[:i], indices[i:]), axis=0)
		indices=permutation(indices)
		a, b, c = indices[0], indices[1], indices[2]
		
		# Get the point (x)
		x=self.population[i]
		
		# Generate mutant (y)
		y=self.population[a]+self.w*(self.population[b]-self.population[c])
		
		# Place it inside the box
		self.bound(y)
		
		# Generate ndim random numbers from 0..1
		uvec=random(self.ndim)
		
		# Crossover x and y, use components of y with probability self.pc
		yind=where(uvec<self.pc)[0]
		z=x.copy()
		z[yind]=y[yind]
		
		return z

def main():
	f = Frame()
	f.pack(side = 'top', expand = 1)
	quit = Button(f, text = 'Quit', command = sys.exit)
	quit.pack(side = 'top')
	o = Opengl(width = 400, height = 400, double = 1)
	a = arange(0,n)
	vertex = shuffle(cos(2*pi*a/n), sin(2*pi*a/n))
	vertex.shape = (n, 2)
#	vertex1 = shuffle(0.5*cos(2*pi*a/n), 0.5*sin(2*pi*a/n))
#	color=ones((n, 3), 'i')
#	color[0]=[1,0,0]
#	color[1]=[1,1,0]
#	color[1]=[1,0,0]
	color = random(n*3)
	color.shape = (n, 3)

	glVertexPointerd(vertex)
	glColorPointerd(color)
	glEnableClientState(GL_VERTEX_ARRAY)
	glEnableClientState(GL_COLOR_ARRAY)

	o.redraw = redraw
	o.pack(side = 'top', expand = 1, fill = 'both')
	o.mainloop()

def find_in_list(distances, sizes, exposure, within_radius):
    """
    Corresponds to the check_and_choose function.
    Given a unit, ask which exposure might accept it as a target.
    """
    for target_idx in within_radius:
        for exp in exposure:
            difference=abs(distances[target_idx]-
                    exp.exposure_distance)
            if exp.target_idx is None:
                exp.target_idx=target_idx
                exp.target_distance=distances[target_idx]
                exp.difference=difference
                exp.cumulative_size+=sizes[target_idx]
            elif abs(difference-exp.difference)<scenario.epsilon:
                chance=np.random()< scenario.size[target_idx]/exp.cumulative_size
                allowed=zoned[target_idx] is False and quarantined[target_idx] is False
                if allowed and chance:
                    exp.target_idx=target_idx
                    exp.target_distance=distances[target_idx]
                    exp.difference=difference
                    exp.cumulative_size+=sizes[target_idx]
            elif difference < exp.difference:
                exp.target_idx=target_idx
                exp.target_distance=distances[target_idx]
                exp.difference=difference
                exp.cumulative_size+=sizes[target_idx]
            else:
                pass

def subsample_imbalanced(X_enhancer, X_promoter, y, positive_subsample_frac):
    n = np.shape(y_train)[0] # sample size (i.e., number of pairs)

    # indices that are positive and selected to be retained or negative
    to_keep = (np.random(n) < positive_subsample_frac) or (y == 1)

    return X_enhancer[to_keep, :], X_promoter[to_keep, :], y[to_keep]

def sim_mh():
    # MH
    x0 = random(), random()
    x = x0
    z = []
    M = 100000
    c = 0

    for i in range(0, M):
        accept, x = mh_step(x)
        if accept:
            c += 1
        z.append(x)

    print("MH Acceptance rate:", 100 * c / M, "%")
    plt.hexbin([x for x, y in z], [y for x, y in z])

def layerModes():
    N = mypaintlib.TILE_SIZE

    dst = np.zeros((N, N, 4), 'uint16')  # rgbu
    dst_values = []
    r1 = range(0, 20)
    r2 = range((1 << 15)/2-10, (1 << 15)/2+10)
    r3 = range((1 << 15)-19, (1 << 15)+1)
    dst_values = r1 + r2 + r3

    src = np.zeros((N, N, 4), 'int64')
    alphas = np.hstack((
        np.arange(N/4),                     # low alpha
        (1 << 15)/2 - np.arange(N/4),       # 50% alpha
        (1 << 15) - np.arange(N/4),         # high alpha
        np.randint((1 << 15)+1, size=N/4),  # random alpha
        ))
    #plot(alphas); show()
    src[:, :, 3] = alphas.reshape(N, 1)  # alpha changes along y axis

    src[:, :, 0] = alphas  # red
    src[:, N*0/4:N*1/4, 0] = np.arange(N/4)  # dark colors
    src[:, N*1/4:N*2/4, 0] = alphas[N*1/4:N*2/4]/2 + np.arange(N/4) - N/2  # 50% lightness
    src[:, N*2/4:N*3/4, 0] = alphas[N*2/4:N*3/4] - np.arange(N/4)  # bright colors
    src[:, N*3/4:N*4/4, 0] = alphas[N*3/4:N*4/4] * np.random(N/4)  # random colors
    # clip away colors that are not possible due to low alpha
    src[:, :, 0] = np.minimum(src[:, :, 0], src[:, :, 3]).clip(0, 1 << 15)
    src = src.astype('uint16')

    #figure(1)
    #imshow(src[:,:,3], interpolation='nearest')
    #colorbar()
    #figure(2)
    #imshow(src[:,:,0], interpolation='nearest')
    #colorbar()
    #show()

    src[:, :, 1] = src[:, :, 0]  # green
    src[:, :, 2] = src[:, :, 0]  # blue

    for name in dir(mypaintlib):
        if not name.startswith('tile_composite_'):
            continue
        junk1, junk2, mode = name.split('_', 2)
        print('testing', name, 'for invalid output')
        f = getattr(mypaintlib, name)
        for dst_value in dst_values:
            for alpha in [1.0, 0.999, 0.99, 0.90, 0.51, 0.50, 0.49, 0.01, 0.001, 0.0]:
                dst[:] = dst_value
                dst_has_alpha = False
                src_opacity = alpha
                f(src, dst, dst_has_alpha, src_opacity)
                #imshow(dst[:,:,0], interpolation='nearest')
                #gray()
                #colorbar()
                #show()
                errors = dst > (1 << 15)
                assert not errors.any()
        print('passed')