def draw(self):
        """
        if sample_priors = True and random_sample = True:
           draw returns a random draw of a categorical distribution with parameters drawn from a Dirichlet distribution
           the hyperparameters on the Dirichlet are given by the bandit's metric with laplacian smoothing
        if sample_priors = False and random_sample = True:
            draw returns a random draw of a categorical distribution with parameters given by the bandit's metric
        if sample_priors = True and random_sample = False:
            draw returns argmax(random.dirichlet((x_0 + 1, ... , x_n_arms + 1))) where x_i is the ith value returned by
            the bandit's metric.
        if sample_priors = False and random_sample = False:
            become a purely greedy bandit with the selected arm given by argmax(metric)

        :return: The numerical index of the selected arm
        """
        temp = self._schedule_fn(self.total_draws)
        x = array(self._metric_fn()) * temp + 1

        if self.sample_priors:
            pvals = random.dirichlet(x)
        else:
            pvals = x / sum(x)

        if self.random_sample:
            return argmax(random.multinomial(1, pvals=pvals))
        else:
            return argmax(pvals)

	def testPercentErrorIsSame(self):
		NN.pat = zip(self.trn_d['input'], self.trn_d['target'])		
		pyb_ws = self.net.params.copy()
		nn = NN()
		nn.wi = pyb_ws[:nn.wi.size].reshape(NN.nh, NN.ni).T
		nn.wo = pyb_ws[nn.wi.size:].reshape(NN.no, NN.nh).T
		correct = 0
		wrong = 0
		argmax_cor = 0
		argmax_wng = 0
		all_aos = []
		for i, x in enumerate(self.trn_d['input']):
			nn.activate(x)
			out = self.net.activate(x)
			# print 'ga bp trg', nn.ao, out, self.trn_d['target'][i], '++++' if not (out - self.trn_d['target'][i]).any() else '-'
			all_aos.append(nn.ao.copy())
			if not (out - self.trn_d['target'][i]).any():
				correct += 1
			else:
				wrong += 1
			if argmax(out) == argmax(self.trn_d['target'][i]):
				argmax_cor += 1
			else:
				argmax_wng += 1
		print 'actual', wrong, 'wrong', correct, 'correct', float(wrong) / (wrong + correct) * 100
		print 'using argmax', argmax_wng, 'wrong', argmax_cor, 'correct', float(argmax_wng) / (argmax_wng + argmax_cor) * 100
		argmax_perc_err = float(argmax_wng) / (argmax_wng + argmax_cor) * 100
		res = nn.sumErrors()
		nn_perc_err = 100 - res[1]
		pb_nn_perc_err = percentError(self.trainer.testOnClassData(), self.trn_d['class'])
		self.assertAlmostEqual(nn_perc_err, pb_nn_perc_err)
		self.assertAlmostEqual(nn_perc_err, pb_nn_perc_err, argmax_perc_err)

def plot_disc_policy():
    #First compute policy function...==========================================
    N = 500
    w = sp.linspace(0,100,N)
    w = w.reshape(N,1)
    u = lambda c: sp.sqrt(c)
    util_vec = u(w)
    alpha = 0.5
    alpha_util = u(alpha*w)
    alpha_util_grid = sp.repeat(alpha_util,N,1)
    
    m = 20
    v = 200
    f = discretelognorm(w,m,v)
    
    VEprime = sp.zeros((N,1))
    VUprime    = sp.zeros((N,N))
    EVUprime = sp.zeros((N,1))
    psiprime = sp.ones((N,1))
    gamma = 0.1
    beta = 0.9
    
    m = 15
    tol = 10**-9
    delta = 1+tol
    it = 0
    while (delta >= tol):
        it += 1
        
        psi = psiprime.copy()
        arg1 = sp.repeat(sp.transpose(VEprime),N,0)
        arg2 = sp.repeat(EVUprime,N,1)
        arg = sp.array([arg2,arg1])
        psiprime = sp.argmax(arg,axis = 0) 
        
        for j in sp.arange(0,m):
            VE = VEprime.copy()
            VU = VUprime.copy()
            EVU = EVUprime.copy()
            VEprime = util_vec + beta*((1-gamma)*VE + gamma*EVU)
            arg1 = sp.repeat(sp.transpose(VE),N,0)*psiprime
            arg2 = sp.repeat(EVU,N,1)*(1-psiprime)
            arg = arg1+arg2
            VUprime = alpha_util_grid + beta*arg
            EVUprime = sp.dot(VUprime,f)  
    
        
    
        delta = sp.linalg.norm(psiprime -psi) 

    wr_ind = sp.argmax(sp.diff(psiprime), axis = 1)
    wr = w[wr_ind]
    print w[250],wr[250]
        
    #Then plot=================================================================
    plt.plot(w,psiprime[250,:]) 
    plt.ylim([-.5,1.5])      
    plt.xlabel(r'$w\prime$')
    plt.yticks([0,1])
    plt.savefig('disc_policy.pdf')

def bandpower(f, Pxx, fmin, fmax):
    """ integrate the power spectral density between fmin and fmax
        using the trapezoidal method
    """
    ind_min = scipy.argmax(f > fmin) - 1
    ind_max = scipy.argmax(f > fmax) - 1
    return scipy.trapz(Pxx[ind_min: ind_max], f[ind_min: ind_max])

 def plot_REL_ERR_SU2(self,which_case):
      i=0;
      thermo1 = self.select[which_case][0]
      thermo2 = self.select[which_case][1]
      get_REL_ERR_SU2(self,which_case)
      
      print 'Median error SU2', sp.median(self.REL_ERR)
      print 'Mean error SU2', sp.mean(self.REL_ERR)
      print 'Max error SU2', max(self.REL_ERR)
      print 'Min error SU2', min(self.REL_ERR)
      x = getattr(self.SU2[which_case],thermo1)
      y = getattr(self.SU2[which_case],thermo2)
      #trusted_values = sp.where(self.REL_ERR>0<0.9*max(self.REL_ERR))
      self.REL_ERR = self.REL_ERR[trusted_values]
      x = x[trusted_values]
      y = y[trusted_values]
      scat=plt.scatter(x,y,c=self.REL_ERR, s=1)                
      plt.grid(which='both')
      scat.set_array(self.REL_ERR)        
      plt.colorbar(scat)
      plt.xlim((min(x)*0.95,max(x)*1.05));
      plt.ylim((min(y)*0.95,max(y)*1.05));
      print 'x argmax %i , x_val: %f ' %(sp.argmax(self.REL_ERR),x[sp.argmax(self.REL_ERR)])
      print 'y argmax %i , y_val: %f ' %(sp.argmax(self.REL_ERR),y[sp.argmax(self.REL_ERR)])
      return;

    def error(self, results, expected):

        err = 0
        for i in range(results.shape[1]):

            err += self.lsexp(results[:, i]) - sp.dot(expected[:, i], results[:, i])

        misclassified = sp.sum(sp.argmax(results, axis=0) != sp.argmax(expected, axis=0))

        return err, misclassified

def Problem6Real():
    N = 500
    w = sp.linspace(0,100,N)
    w = w.reshape(N,1)
    u = lambda c: sp.sqrt(c)
    util_vec = u(w)
    alpha = 0.5
    alpha_util = u(alpha*w)
    alpha_util_grid = sp.repeat(alpha_util,N,1)
    
    m = 20
    v = 200
    f = discretelognorm(w,m,v)
    
    VEprime = sp.zeros((N,1))
    VUprime    = sp.zeros((N,N))
    EVUprime = sp.zeros((N,1))
    psiprime = sp.ones((N,1))
    gamma = 0.1
    beta = 0.9
    
    m = 15
    tol = 10**-9
    delta = 1+tol
    it = 0
    while (delta >= tol):
        it += 1
        
        psi = psiprime.copy()
        arg1 = sp.repeat(sp.transpose(VEprime),N,0)
        arg2 = sp.repeat(EVUprime,N,1)
        arg = sp.array([arg2,arg1])
        psiprime = sp.argmax(arg,axis = 0)    
        
        for j in sp.arange(0,m):
            VE = VEprime.copy()
            VU = VUprime.copy()
            EVU = EVUprime.copy()
            VEprime = util_vec + beta*((1-gamma)*VE + gamma*EVU)
            arg1 = sp.repeat(sp.transpose(VE),N,0)*psiprime
            arg2 = sp.repeat(EVU,N,1)*(1-psiprime)
            arg = arg1+arg2
            VUprime = alpha_util_grid + beta*arg
            EVUprime = sp.dot(VUprime,f)  
    
        
    
        delta = sp.linalg.norm(psiprime -psi)
        #print(delta)    
        
    wr_ind = sp.argmax(sp.diff(psiprime), axis = 1)
    wr = w[wr_ind]
    plt.plot(w,wr)
    plt.show()
    return wr

def generate(hmm,observation_space,n_sim):
	A = hmm[0]
	B = hmm[1]
	pi = hmm[2]
	states = sp.zeros(n_sim)
	observations = []
	states[0] = sp.argmax(sp.random.multinomial(1,pi))
	observations.append(observation_space[sp.argmax(sp.random.multinomial(1,B[states[0],:]))])
	for i in range(1,n_sim):
		states[i] = sp.argmax(sp.random.multinomial(1,A[states[i-1],:]))
		observations.append(observation_space[sp.argmax(sp.random.multinomial(1,B[states[i],:]))])
	return states,observations

def pitch(x, fs, pitchrange=[12,120], mode='corr'):
    if mode=='corr':
        corr = scipy.correlate(x, x, mode='full')[len(x)-1:]
        corr[:int(fs/midi2hz(pitchrange[1]))] = 0
        corr[int(fs/midi2hz(pitchrange[0])):] = 0
        indmax = scipy.argmax(corr)
    elif mode=='ceps':
        y = rceps(x)
        y[:int(fs/midi2hz(pitchrange[1]))] = 0
        y[int(fs/midi2hz(pitchrange[0])):] = 0
        indmax = scipy.argmax(y)
    return hz2midi(fs/indmax)

def calcError(trainer, dataset=None):
    if dataset == None:
        dataset = trainer.ds
    dataset.reset()
    out = []
    targ = []
    for seq in dataset._provideSequences():
        trainer.module.reset()
        for input, target in seq:
            res = trainer.module.activate(input)
            out.append(argmax(res))
            targ.append(argmax(target))
    return percentError(out, targ) / 100

def generateGaussianHMM(hmm,n_sim):
	A = hmm[0]
	means = hmm[1]
	covars = hmm[2]
	pi = hmm[3]
	states = sp.zeros(n_sim).astype(int)
	K = len(means[0,:])
	observations = sp.zeros((n_sim,K))
	states[0] = int(sp.argmax(sp.random.multinomial(1,pi)))
	observations[0,:] = sp.random.multivariate_normal(means[states[0],:],covars[states[0],:,:])
	for i in range(1,n_sim):
		states[i] = int(sp.argmax(sp.random.multinomial(1,A[states[i-1],:])))
		observations[i,:] = sp.random.multivariate_normal(means[states[i],:],covars[states[i],:,:])
	return states,observations

def digshc(docs, alpha, threshold, epsilon, hr_min):
    # Worth nothing that this takes in plaintext documents, _not_ vectors.
    shc = SimilarityHistogramClusterer(alpha, threshold, epsilon, hr_min)

    for doc in docs:
        shc.fit(doc)

    doc_clus_map = {}
    for idx, clus in enumerate(shc.formed_clusters):
        for doc_id in clus.doc_ids:
            doc_clus_map.setdefault(doc_id, [])
            doc_clus_map[doc_id].append(idx)

    labels = []
    for id in sorted(doc_clus_map):
        cluster_ids = doc_clus_map[id]
        if len(cluster_ids) == 1:
            best_cl_id = cluster_ids[0]
        else:
            clusters = [shc.get_cluster(cl_id) for cl_id in cluster_ids]
            sims = [shc.get_cluster_sim(cl, shc.get_doc(id)) for cl in clusters]
            max_i = argmax(sims)
            best_cl_id = clusters[max_i].id
        labels.append(best_cl_id)

    return labels

 def testKernelCoeffs(self):
     for scale in [0.35, 0.5, 0.75, 1]:
         for dim in [0,1,2,3,4]:
             dgK = mango.image.discrete_gaussian_kernel(sigma=scale, dim=dim, errtol=0.001)
             self.assertAlmostEqual(1.0, sp.sum(dgK), 8)
             mxElem = sp.argmax(dgK)
             self.assertTrue(sp.all(sp.array(dgK.shape)//2 == sp.unravel_index(mxElem, dgK.shape)))

def getPredominantColor(filename):
    im = Image.open(filename).convert('RGB')

    # Convert to numpy array
    ar = scipy.misc.fromimage(im)

    # Get dimensions
    shape = ar.shape

    # Convert to bidimensional array of width x height rows and 3 columns (RGB)
    ar = ar.reshape(scipy.product(shape[:2]), shape[2])

    # Find cluster centers and their distortions
    # codes contains the RGB value of the centers
    codes, dist = scipy.cluster.vq.kmeans(ar.astype(float), NUM_CLUSTERS)

    # Maps all the pixels in the image to their respective centers
    vecs, dist = scipy.cluster.vq.vq(ar, codes)

    # Counts the occurances of each color (NUM_CLUSTER different colors after the mapping)
    counts, bins = scipy.histogram(vecs, len(codes))

    # Find most frequent color
    index_max = scipy.argmax(counts)
    peak = codes[index_max]

    return peak.astype(int)

def getDominantColor(img_url):
    if r.exists(img_url):
        cache_result = r.hmget(img_url, ['r', 'g', 'b'])
        return cache_result
        
    NUM_CLUSTERS = 5
    im = Image.open(StringIO.StringIO(urllib2.urlopen(img_url).read()))
    img_arr = scipy.misc.fromimage(im)
    img_shape = img_arr.shape
    
    if len(img_shape) > 2:
        img_arr = img_arr.reshape(scipy.product(img_shape[:2]), img_shape[2])
    
    codes, _ = scipy.cluster.vq.kmeans(img_arr, NUM_CLUSTERS)
    
    original_codes = codes
    for low, hi in [(60, 200), (35, 230), (10, 250)]:
        codes = scipy.array([code for code in codes if not (all([c < low for c in code]) or all([c > hi for c in code]))])
        if not len(codes):
            codes = original_codes
        else:
            break

    vecs, _ = scipy.cluster.vq.vq(img_arr, codes)
    counts, bins = scipy.histogram(vecs, len(codes))

    index_max = scipy.argmax(counts)
    peak = codes[index_max]
    color = [int(c) for c in peak[:3]]
    r.hmset(img_url, {'r':color[0], 'g':color[1], 'b':color[2]})
    #r.expire(img_url, 86400)
    return color

def displayData(X, theta = None):
    """Display 2D data in a nice grid"""
    width = 20
    rows, cols = 10, 10
    out = sp.zeros((width * rows, width * cols))

    rand_indices = sp.random.permutation(5000)[0:rows * cols]

    counter = 0
    for y in range(0, rows):
        for x in range(0, cols):
            start_x = x * width
            start_y = y * width
            out[start_x:start_x+width, start_y:start_y+width] = X[rand_indices[counter]].reshape(width, width).T
            counter += 1

    img = sp.misc.toimage(out)
    figure = plt.figure()
    axes = figure.add_subplot(111)
    axes.imshow(img)

    if theta is not None:
        result_matrix = []
        X_biased = sp.c_[sp.ones(X.shape[0]), X]

        for idx in rand_indices:
            result = (sp.argmax(theta.T.dot(X_biased[idx])) + 1) % 10
            result_matrix.append(result)
        result_matrix = sp.array(result_matrix).reshape(rows, cols).transpose()
        print(result_matrix)

    plt.show()

def Problem3Real():
    beta = 0.9
    N = 1000
    u = lambda c: sp.sqrt(c)
    W = sp.linspace(0,1,N)
    X, Y = sp.meshgrid(W,W)
    Wdiff = sp.transpose(X-Y)
    index = Wdiff <0
    Wdiff[index] = 0
    util_grid = u(Wdiff)
    util_grid[index] = -10**10
    
    Vprime = sp.zeros((N,1))
    psi = sp.zeros((N,1))
    delta = 1.0
    tol = 10**-9
    it = 0
    max_iter = 500
    
    while (delta >= tol) and (it < max_iter):
        V = Vprime
        it += 1;
        #print(it)
        val = util_grid + beta*sp.transpose(V)
        Vprime = sp.amax(val, axis = 1)
        Vprime = Vprime.reshape((N,1))
        psi_ind = sp.argmax(val,axis = 1)
        psi    = W[psi_ind]
        delta = sp.dot(sp.transpose(Vprime - V),Vprime-V)
    
    return psi

 def _box_cox_transform(self, verbose=False, method='standard'):
     """
     Performs the Box-Cox transformation, over different ranges, picking the optimal one w. respect to normality.
     """
     from scipy import stats
     a = sp.array(self.values)
     if method == 'standard':
         vals = (a - min(a)) + 0.1 * sp.var(a)
     else:
         vals = a
     sw_pvals = []
     lambdas = sp.arange(-2.0, 2.1, 0.1)
     for l in lambdas:
         if l == 0:
             vs = sp.log(vals)
         else:
             vs = ((vals ** l) - 1) / l
         r = stats.shapiro(vs)
         if sp.isfinite(r[0]):
             pval = r[1]
         else:
             pval = 0.0
         sw_pvals.append(pval)
     i = sp.argmax(sw_pvals)
     l = lambdas[i]
     if l == 0:
         vs = sp.log(vals)
     else:
         vs = ((vals ** l) - 1) / l
     self._perform_transform(vs,"box_cox")
     log.debug('optimal lambda was %0.1f' % l)
     return True

 def generate_threshold_mesh(self, min_value=0.0, max_value=1.0e9):
     r"""
     Generates a mesh excluding all blocks below the min_value arg. Regions
     that are isolated by the thresholding are also automatically removed.
     """
     #
     # thresholding the data and then checking for isolated clusters
     self._field.threshold_data(min_value, max_value, repl=0.0)
     self._field.copy_data(self)
     #
     adj_matrix = self._field.create_adjacency_matrix()
     num_cs, cs_ids = csgraph.connected_components(csgraph=adj_matrix,
                                                   directed=False)
     # only saving the largest cluster
     if num_cs > 1:
         cs_count = sp.zeros(num_cs, dtype=int)
         for cs_num in cs_ids:
             cs_count[cs_num] += 1
         self.data_vector[sp.where(cs_ids != sp.argmax(cs_count))[0]] = 0.0
         self.data_map = sp.reshape(self.data_vector, (self.nz, self.nx))
     #
     self._field.data_map = self.data_map
     self._field.data_vector = sp.ravel(self.data_map)
     #
     # generating blocks and vertices
     mask = self.data_map > 0.0
     self._generate_masked_mesh(cell_mask=mask)

	def classify(self, xL, xR):

		x = sp.vstack([xL,xR, sp.ones(xR.shape[1])]).T

		tmp = sp.dot(x, self.w)

		return tmp, sp.argmax(tmp, axis=1)