def MaxImpedanceComputation(InputGraph):
	
	MaxTotalImpedance=0
	
	G=InputGraph.copy()
	number_of_vertices=G.order()
	vertexlist=G.nodes()
	
	for top_node in vertexlist:
		for ground_node in vertexlist:
			if ground_node<top_node:
				ordered_vertexlist=vertexlist[:]
				ordered_vertexlist.remove(top_node)
				ordered_vertexlist.remove(ground_node)
				ordered_vertexlist.insert(0,top_node)
				ordered_vertexlist.insert(0,ground_node)
				
				LaplacianMatrix=nx.laplacian(G,ordered_vertexlist)
				ConductanceMatrix=np.delete(LaplacianMatrix,0,0)
				ConductanceMatrix=np.delete(ConductanceMatrix,np.s_[0],1)
				InputVector=[0]*(number_of_vertices-1)
				InputVector[0]=1
				VoltageVector=linalg.solve(ConductanceMatrix,InputVector)
				TotalImpedance=VoltageVector[0]
				if TotalImpedance>MaxTotalImpedance:
					MaxTotalImpedance=TotalImpedance
	
	return MaxTotalImpedance

def non_max_suppression(boxes, scores, threshold):
    """Performs non-maximum supression and returns indicies of kept boxes.
    boxes: [N, (y1, x1, y2, x2)]. Notice that (y2, x2) lays outside the box.
    scores: 1-D array of box scores.
    threshold: Float. IoU threshold to use for filtering.
    """
    assert boxes.shape[0] > 0
    if boxes.dtype.kind != "f":
        boxes = boxes.astype(np.float32)

    # Compute box areas
    y1 = boxes[:, 0]
    x1 = boxes[:, 1]
    y2 = boxes[:, 2]
    x2 = boxes[:, 3]
    area = (y2 - y1) * (x2 - x1)

    # Get indicies of boxes sorted by scores (highest first)
    ixs = scores.argsort()[::-1]

    pick = []
    while len(ixs) > 0:
        # Pick top box and add its index to the list
        i = ixs[0]
        pick.append(i)
        # Compute IoU of the picked box with the rest
        iou = compute_iou(boxes[i], boxes[ixs[1:]], area[i], area[ixs[1:]])
        # Identify boxes with IoU over the threshold. This
        # returns indicies into ixs[1:], so add 1 to get
        # indicies into ixs.
        remove_ixs = np.where(iou > threshold)[0] + 1
        # Remove indicies of the picked and overlapped boxes.
        ixs = np.delete(ixs, remove_ixs)
        ixs = np.delete(ixs, 0)
    return np.array(pick, dtype=np.int32)

    def append_new_point(self, y, x=None):
        self._axis_y_array = np.append(self._axis_y_array, y)
        if x:
            self._axis_x_array = np.append(self._axis_x_array, x)
        else:
            self._axis_x_array = np.arange(len(self._axis_y_array))

        if self.max_plot_points:
            if self._axis_y_array.size > self.max_plot_points:
                self._axis_y_array = np.delete(self._axis_y_array, 0)
                self._axis_x_array = np.delete(self._axis_x_array, 0)

        if self.single_curve is None:
            self.single_curve, = self.axes.plot(
                self._axis_y_array, linewidth=2, marker="s"
            )
        else:
            self.axes.fill(self._axis_y_array, "r", linewidth=2)

        self._axis_y_limits[1] = (
            self._axis_y_array.max() + self._axis_y_array.max() * 0.05
        )
        self.axes.set_ylim(self._axis_y_limits)
        self.single_curve.set_xdata(self._axis_x_array)
        self.single_curve.set_ydata(self._axis_y_array)
        self.axes.relim()
        self.axes.autoscale_view()
        self.fig.canvas.draw()
        self.fig.canvas.flush_events()
        self.axes.grid(True)

        # TODO move y lims as propery
        self.axes.set_ylim(
            (0, self._axis_y_array.max() + self._axis_y_array.max() * 0.05)
        )

def carbonylorcarboxyl(allligand,index,bond_dist):

	allligandcoods = allligand.positions
	ocoods = np.zeros((1,3), dtype = float)
	ocoods[0,:] = allligandcoods[index,:]
	ocoods = np.float32(ocoods)

	tempdist = MDAnalysis.lib.distances.distance_array(ocoods,allligandcoods)
	A = np.argsort(tempdist)
	temp = int(A[0,1])

	Omatecood = np.zeros((1,3), dtype = float)
	Omatecood[0,:] = allligandcoods[temp,:]
	Omatecood = np.float32(Omatecood)

	tempdist2 = MDAnalysis.lib.distances.distance_array(Omatecood, allligandcoods)
	B = np.argsort(tempdist2)
	B = np.delete(B,0,axis = 1)
	for i in xrange(0,B.size):
		if B[0,i] == index:
			C = np.delete(B,i,axis = 1)
			break

	base1 = int(C[0,0])
	base2 = int(C[0,1])
	type1 = allligand[base1].type
	type2 = allligand[base2].type

	if type1 == 'O' or type2 == 'O':
		atype = 'carboxyl'
	else:
		atype = 'carbonyl'

	return atype

 def make_DeviationPlot(self,year):
     average = np.array(self.average)
     deviation = np.array(self.deviation)
     
     dates = []
     dis = []
     count = 0
     for d in self.time:
         if self.time[count].year == year:
             dates.append(datetime.date(self.time[count].year, self.time[count].month,self.time[count].day))
             dis.append(self.discharge[count])
             
         count += 1
     
     dis = np.array(dis)
     dates = np.array(dates)  
     
     if len(dates) == 365:
         average = np.delete(average,-1)
         deviation = np.delete(deviation,-1)
         
     plus1 = np.array(average+deviation)
     minus1 = np.array(average-deviation)
     
             
     plt.plot(dates,dis,'r')
     x = np.linspace(1,366,366)
     plt.plot(dates,average,'k')
     plt.fill_between(dates,average,plus1,facecolor='gray')
     plt.fill_between(dates,average,minus1,facecolor='gray')

	def StoreTransition(self, s_t, a_t, r_t, s_t_next, d_t=0):
		s_t = s_t.reshape(1, self.state_size)
		s_t_next = s_t_next.reshape(1, self.state_size)
		a_t = a_t.reshape(1, self.action_size)
		r_t = r_t.reshape(1, 1)
		d_t = np.array([d_t]).reshape(1, 1)
		
		self.S = np.concatenate((self.S, s_t))
		self.Stag = np.concatenate((self.Stag, s_t_next))
		self.A = np.concatenate((self.A, a_t))
		self.R = np.concatenate((self.R, r_t))
		self.D = np.concatenate((self.D, d_t))

		if self.populated < self.buffer_size:
			if self.populated == 0:
				self.S = np.delete(self.S,0,0)
				self.A = np.delete(self.A,0,0)
				self.R = np.delete(self.R,0,0)
				self.Stag = np.delete(self.Stag,0,0)
				self.D = np.delete(self.D,0,0)
			self.populated += 1
		else:
			self.S = np.delete(self.S,0,0)
			self.A = np.delete(self.A,0,0)
			self.R = np.delete(self.R,0,0)
			self.Stag = np.delete(self.Stag,0,0)
			self.D = np.delete(self.D,0,0)

 def execEnd(self,eventIdx):
     # execute an end-breaking or depolymerization event.
     oligoEndBreak=self.ald['end'][eventIdx/2]
     leftRight=eventIdx%2*2-1
     lr=-(leftRight+1)/2
     unitMoving=oligoEndBreak.ends[lr]
     oligo_vanish,form_oligo,self.event_code=oligoEndBreak.end_break(leftRight,self.units)
     if form_oligo:
         # not empty
         mono=form_oligo['monomer']
         if mono:
             # monomer + monomer (mergeOligo)
             idx=np.where([x in [mono,unitMoving] for x in self.monomers])[0]
             self.monomers=np.delete(self.monomers,idx)
             self.oligos=np.insert(self.oligos,0,form_oligo['oligo'])
         else:
             # monomer + multimer (mergeOligo)
             idx=np.where([unitMoving is x for x in self.monomers])[0]
             self.monomers=np.delete(self.monomers,idx)
     else:
         #empty, add the end to monomers
         self.monomers=np.insert(self.monomers,0,unitMoving)
         unitMoving.energize()
     
     if oligo_vanish:
         idx=np.where([oligoEndBreak is x for x in self.oligos])[0]
         self.oligos=np.delete(self.oligos,idx)
         
         idx=np.where([unitMoving is not x for x in oligoEndBreak.subunits])[0]
         nonmoving_unit=oligoEndBreak.subunits[idx[0]]
         self.monomers=np.insert(self.monomers,0,nonmoving_unit)
         nonmoving_unit.energize()

def sineModelSynth(tfreq, tmag, tphase, N, H, fs):
	"""
	Synthesis of a sound using the sinusoidal model
	tfreq,tmag,tphase: frequencies, magnitudes and phases of sinusoids
	N: synthesis FFT size, H: hop size, fs: sampling rate
	returns y: output array sound
	"""
	
	hN = N/2                                                # half of FFT size for synthesis
	L = tfreq.shape[0]                                      # number of frames
	pout = 0                                                # initialize output sound pointer         
	ysize = H*(L+3)                                         # output sound size
	y = np.zeros(ysize)                                     # initialize output array
	sw = np.zeros(N)                                        # initialize synthesis window
	ow = triang(2*H)                                        # triangular window
	sw[hN-H:hN+H] = ow                                      # add triangular window
	bh = blackmanharris(N)                                  # blackmanharris window
	bh = bh / sum(bh)                                       # normalized blackmanharris window
	sw[hN-H:hN+H] = sw[hN-H:hN+H]/bh[hN-H:hN+H]             # normalized synthesis window
	lastytfreq = tfreq[0,:]                                 # initialize synthesis frequencies
	ytphase = 2*np.pi*np.random.rand(tfreq[0,:].size)       # initialize synthesis phases 
	for l in range(L):                                      # iterate over all frames
		if (tphase.size > 0):                                 # if no phases generate them
			ytphase = tphase[l,:] 
		else:
			ytphase += (np.pi*(lastytfreq+tfreq[l,:])/fs)*H     # propagate phases
		Y = UF.genSpecSines(tfreq[l,:], tmag[l,:], ytphase, N, fs)  # generate sines in the spectrum         
		lastytfreq = tfreq[l,:]                               # save frequency for phase propagation
		ytphase = ytphase % (2*np.pi)                         # make phase inside 2*pi
		yw = np.real(fftshift(ifft(Y)))                       # compute inverse FFT
		y[pout:pout+N] += sw*yw                               # overlap-add and apply a synthesis window
		pout += H                                             # advance sound pointer
	y = np.delete(y, range(hN))                             # delete half of first window
	y = np.delete(y, range(y.size-hN, y.size))              # delete half of the last window 
	return y

def project_into_plane(index, r0, rm):
    r'''Projects out-of-plane resolution into a specified plane by performing
    a gaussian integral over the third axis.

    Parameters
    ----------
    index : int
        Index of the axis that should be integrated out

    r0 : float
        Resolution prefactor

    rm : ndarray
        Resolution array

    Returns
    -------
    mp : ndarray
        Resolution matrix in a specified plane

    '''

    r = np.sqrt(2 * np.pi / rm[index, index]) * r0
    mp = rm

    b = rm[:, index] + rm[index, :].T
    b = np.delete(b, index, 0)

    mp = np.delete(mp, index, 0)
    mp = np.delete(mp, index, 1)

    mp -= 1 / (4. * rm[index, index]) * np.outer(b, b.T)

    return [r, mp]

def edit_description(instance):

    # twenty different categories
    scores = [0] * 20

    # Strip out all the punctuation
    unstripped = instance[9].lower()
    for c in string.punctuation:
        unstripped = unstripped.replace(c,"")

    description = unstripped.split()

    # add to the score if a word matches a category
    # 10 is the description
    for word in description:
        for i, category in enumerate(LDA):
            if word in category:
                scores[i] += 1

    # save the target
    target = instance[-1]

    # get rid of the description and target columns
    instance = np.delete(instance, 10, 0) # 10 is which column, 1 means column, 0 means row
    instance = np.delete(instance, -1, 0)

    # add the scores
    instance = np.append(instance, scores)

    # add the target back on the end
    return np.append(instance, target)

def build_tree(data, labels, word_data, level):
    if (level == 0):
        #return label value which is dominant
        return LabelConv[st.mode(labels)[0][0]-1];
    #select appropriate attribute for the node:
    best, best_ig = attribute_selection(data,labels);
    best_data = data[:,best]; best_word = word_data[best];
    #remove all regarding that attribute from the data:
    word_data = np.delete(word_data,best,0);
    left_data = np.delete(data[best_data == 0,:],best,1); 
    right_data = np.delete(data[best_data == 1,:],best,1);
    #divide labels into two subarray based on selected attribute:
    left_labl = labels[best_data == 0]; 
    right_labl = labels[best_data == 1];
    if (check_label(left_labl) == 2 and level != 0):
        #since label is mono-valued:
        left = LabelConv[left_labl[0]-1];
    else:
        left = build_tree(left_data,left_labl,word_data,level-1);
    if (check_label(right_labl) == 2 and level != 0):
        #since label is mono-valued:
        right = LabelConv[right_labl[0]-1];
    else:
        right = build_tree(right_data,right_labl,word_data,level-1);
    subtrees = {0: left, 1: right};
    return (best_word,best_ig,subtrees);

def pixel_coordinates(nx, ny, mode="centers"):
    """Get pixel coordinates from a regular grid with dimension nx by ny.

    Parameters
    ----------
    nx : int
        xsize
    ny : int
        ysize
    mode : string
        `centers` or `edges` to return the pixel coordinates
        defaults to centers
    Returns
    -------
    coordinates : :class:`numpy:numpy.ndarray`
         Array of shape (ny,nx) with pixel coordinates (x,y)

    """
    x = np.linspace(0, nx, num=nx + 1)
    y = np.linspace(0, ny, num=ny + 1)
    if mode == "centers":
        x = x + 0.5
        y = y + 0.5
        x = np.delete(x, -1)
        y = np.delete(y, -1)
    X, Y = np.meshgrid(x, y)
    coordinates = np.empty(X.shape + (2,))
    coordinates[:, :, 0] = X
    coordinates[:, :, 1] = Y
    return coordinates

def reduce_dimension(m):
    """
    reduce the dimension of game matrix based on domination --
    player I will be better off if one row constently larger than another
    player II will be better off if one col constently smaller than anthoer
    Output: the reduced-size game matrix
    Note: This implements stric domination.
    TODO: convex reduction
    """
    local = np.array(m)
    flag = True
    while True:
        rbefore = len(local)
        candidates = []
        for nr in permutations(range(len(local)), 2):
            bigger = reduce(lambda x,y: x and y, local[nr[0]]>local[nr[1]])
            if bigger: 
                candidates.append(nr[1])
        for i in candidates:
            local = np.delete(local, i, 0)

        cbefore = len(local[0])
        candidates = []
        for nc in permutations(range(len(local[0])), 2):
            smaller = reduce(lambda x,y: x and y, local[:,nc[0]]<local[:, nc[1]])
            if smaller:
                candidates.append(nc[1])
        for i in candidates:
            local = np.delete(local, i, 1)

        if len(local[0])==cbefore and len(local)==rbefore:
            break

    return local

def stochasticModelSynth(stocEnv, H, N):
	"""
	Stochastic synthesis of a sound
	stocEnv: stochastic envelope; H: hop size; N: fft size
	returns y: output sound
	"""

	if not(UF.isPower2(N)):                                 	# raise error if N not a power of two
		raise ValueError("N is not a power of two")
 
	hN = N/2+1                                            		# positive size of fft
	No2 = N/2							# half of N
	L = stocEnv[:,0].size                                    	# number of frames
	ysize = H*(L+3)                                         	# output sound size
	y = np.zeros(ysize)                                     	# initialize output array
	ws = 2*hanning(N)                                        	# synthesis window
	pout = 0                                                 	# output sound pointer
	for l in range(L):                    
		mY = resample(stocEnv[l,:], hN)                        # interpolate to original size
		pY = 2*np.pi*np.random.rand(hN)                        # generate phase random values
		Y = np.zeros(N, dtype = complex)                       # initialize synthesis spectrum
		Y[:hN] = 10**(mY/20) * np.exp(1j*pY)                   # generate positive freq.
		Y[hN:] = 10**(mY[-2:0:-1]/20) * np.exp(-1j*pY[-2:0:-1]) # generate negative freq.
		fftbuffer = np.real(ifft(Y))                           # inverse FFT
		y[pout:pout+N] += ws*fftbuffer                         # overlap-add
		pout += H  
	y = np.delete(y, range(No2))                              # delete half of first window
	y = np.delete(y, range(y.size-No2, y.size))               # delete half of the last window 
	return y

  def update_proximity_matrix(self, old_prox, new_centroid, a, b):
    old_prox = np.delete(old_prox, [a,b], 0) #delete rows
    old_prox = np.delete(old_prox, [a,b], 1) #delete cols

    # add a line of zeroes on the right and bottom edges
    mid = np.hstack((old_prox, np.zeros((old_prox.shape[0], 1), dtype=old_prox.dtype)))
    pprint(("mid", mid, mid.shape))
    new_prox = np.vstack((mid, np.zeros((1, mid.shape[1]), dtype=mid.dtype)))

    pprint(("expanded", new_prox, new_prox.shape))

    old_length = len(old_prox) - 1
    new_length = len(new_prox) - 1

    #fill them in with new comparisons
    new_prox[new_length,new_length] = float(HIGH)

    for i, centroid in enumerate(self.centroids[:-1]):
      diff = np.linalg.norm(centroid - new_centroid)

      pprint(("checking", diff, i))

      new_prox[new_length,i] = diff
      new_prox[i,new_length] = diff

    pprint(("new prox", new_prox, new_prox.shape))

    return new_prox

def update_extra_mat(matfile,to_remove):
    """ updates the time_frames, confounds and mask_suppressed arrays to
    reflect the removed volumes. However, does not change other items in
    _extra.mat file

    """

    mat = loadmat(matfile)
    # update time_frames
    ntf = np.delete(mat['time_frames'][0],to_remove)
    mat.update({'time_frames': ntf})

    # update confounds
    ncon = np.delete(mat['confounds'],to_remove,axis = 0)
    mat.update({'confounds': ncon})

    # update mask_suppressed
    ms = mat['mask_suppressed']
    for supp in to_remove:
        ms[supp][0] = 1
    mat.update({'mask_suppressed': ms})

    # save updated mat file
    jnk, flnme = os.path.split(matfile)
    savemat(os.path.join(output_dir,flnme),mat)

    def _rebuild_iso(self, sel):
        g = self.graph
        ss = [p.plots[pp][0] for p in g.plots
              for pp in p.plots
              if pp == 'data{}'.format(self.group_id)]

        self._set_renderer_selection(ss, sel)

        if self._cached_data:
            reg=self._cached_reg
            xs, ys, xerr, yerr = self._cached_data

            nxs = delete(xs, sel)
            nys = delete(ys, sel)
            nxerr = delete(xerr, sel)
            nyerr = delete(yerr, sel)

            # reg = ReedYorkRegressor(xs=nxs, ys=nys,
            #                         xserr=nxerr, yserr=nyerr)
            reg.trait_set(xs=nxs, ys=nys,xserr=nxerr, yserr=nyerr)
            reg.calculate()

            fit = self.graph.plots[0].plots['fit{}'.format(self.group_id)][0]

            mi, ma = self.graph.get_x_limits()
            rxs = linspace(mi, ma, 10)

            rys = reg.predict(rxs)

            fit.index.set_data(rxs)
            fit.value.set_data(rys)

            if self._plot_label:
                self._add_info(self.graph.plots[0], reg, label=self._plot_label)

def extract_sub_matrix(mat, inds):
    """ Extract submatrix of `mat` by deleting `inds` rows/cols
    """
    for i in sorted(inds, reverse=True):
        mat = np.delete(mat, i, axis=0)
        mat = np.delete(mat, i, axis=1)
    return mat

 def loadGlob(self, simu, Z, S):
     file_root = self._fileRoot(simu, Z, S)
     data = np.loadtxt(self.dir + file_root + '/' + file_root + '_all.deus_histo.txt')
     data = np.delete(data, self.nb_histo-1)        
     densityscale = np.linspace(self.glob_start,self.glob_end,self.nb_histo,0)
     densityscale = np.delete(densityscale, self.nb_histo-1)        
     return densityscale, data

    def model_and_predict(self, X_train, y_train, X_test):
        district_idx = self.columns.index('PdDistrict')
        districts = set(X_train[:,district_idx])
        district_ys = {}
        # Grow forest and predict separately for each district's records
        for d in districts:
            district_X_train = X_train[X_train[:, district_idx] == d]
            district_X_train = np.delete(district_X_train, district_idx, 1)
            district_y_train = y_train[X_train[:, district_idx] == d]
            district_X_test = X_test[X_test[:, district_idx] == d]
            district_X_test = np.delete(district_X_test, district_idx, 1)
            print "Growing forest for", d

            # Not saving output in Git so make this deterministic 
            # with random_state
            rf = RandomForestClassifier(n_estimators=self.n_trees, n_jobs=-1,
                                        random_state=782629)
            rf.fit(district_X_train, district_y_train)

            district_ys[d] = list(rf.predict(district_X_test))
            print "Finished", d

        print "All predictions made"

        y_hat = []
        for row in X_test:
            d_ys = district_ys[row[district_idx]]
            y_hat.append(d_ys.pop(0))

        return y_hat