def calcFreqs(X, timeStep, minFreq=0, maxFreq=np.inf):
    if (X.ndim > 1):
        freqs = scipy.fftpack.fftfreq(X.shape[1], timeStep)
        idx1 = np.argsort(freqs)
        freqs = freqs[idx1]
        idx2 = np.where((freqs >= minFreq) & (freqs <= maxFreq))[0]
        freqs = freqs[idx2]
        return freqs, idx1, idx2, 0
    else:
        # sometimes no all the time steps are the same
        allFreqs, lengths = [], []
        idx1s, idx2s = [], []
        if (isinstance(timeStep, float)):
            timeStep = np.ones((len(X))) * timeStep
        for x, dt in zip(X, timeStep):
            freqs = scipy.fftpack.fftfreq(x.shape[0], dt)
            idx1 = np.argsort(freqs)
            freqs = freqs[idx1]
            idx2 = np.where((freqs > minFreq) & (freqs < maxFreq))[0]
            freqs = freqs[idx2]
            allFreqs.append(freqs)
            lengths.append(len(freqs))
            idx1s.append(idx1)
            idx2s.append(idx2)
        maxLenInd = np.argmax(lengths)
        return allFreqs, idx1s, idx2s, maxLenInd

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 _test_corr(old_func, new_func, sel_item):
    from nose.tools import assert_equal, assert_raises
    n_obs = 20
    n_dims = 10
    np.random.seed(0)
    y = np.random.rand(n_obs) * n_obs
    X = np.tile(y, [n_dims, 1]).T + np.random.randn(n_obs, n_dims)
    rho_fast = new_func(X, y)
    # test dimensionality
    assert_equal(rho_fast.ndim, 1)
    assert_equal(rho_fast.shape[0], n_dims)
    # test data
    rho_slow = np.ones(n_dims)
    for dim in range(n_dims):
        rho_slow[dim] = np.array(old_func(X[:, dim], y)).item(sel_item)
    np.testing.assert_array_equal(rho_fast.shape, rho_slow.shape)
    np.testing.assert_array_almost_equal(rho_fast, rho_slow)
    # test errors
    new_func(np.squeeze(X[:, 0]), y)
    assert_raises(ValueError, new_func, y, X)
    assert_raises(ValueError, new_func, X, y[1:])
    # test dtype
    X = np.argsort(X, axis=0) * 2  # ensure no bug at normalization
    y = np.argsort(y, axis=0) * 2
    rho_fast = new_func(X, y, dtype=int)
    rho_slow = np.ones(n_dims)
    for dim in range(n_dims):
        rho_slow[dim] = np.array(old_func(X[:, dim], y)).item(sel_item)
    np.testing.assert_array_almost_equal(rho_fast, rho_slow)

def rankImages( imdists, query_id, dist_type ):
    # PRE [DO NOT TOUCH]
    ranking = []

    # WRITE YOUR CODE HERE
    related_img = []
    related_img = imdists[query_id,:]
    
    # smaller, order asc
    if dist_type == 'euclidean':
        ranking = np.argsort(related_img)
        
    # larger, order desc
    elif dist_type == 'l2':
        ranking = np.argsort(-related_img)
    
    # larger, order desc
    elif dist_type == 'intersect' or dist_type == 'l1':
        ranking = np.argsort(-related_img)
        
    # smaller, order asc
    elif dist_type == 'chi2':
        ranking = np.argsort(related_img)

    # larger, order desc
    elif dist_type == 'hellinger':
        ranking = np.argsort(-related_img)
    
    
    # RETURN [DO NOT TOUCH]
    return ranking

def trustworthiness(X, X_embedded, n_neighbors=5, precomputed=False):
    """Expresses to what extent the local structure is retained.

    The trustworthiness is within [0, 1]. It is defined as

    .. math::

        T(k) = 1 - \frac{2}{nk (2n - 3k - 1)} \sum^n_{i=1}
            \sum_{j \in U^{(k)}_i} (r(i, j) - k)

    where :math:`r(i, j)` is the rank of the embedded datapoint j
    according to the pairwise distances between the embedded datapoints,
    :math:`U^{(k)}_i` is the set of points that are in the k nearest
    neighbors in the embedded space but not in the original space.

    * "Neighborhood Preservation in Nonlinear Projection Methods: An
      Experimental Study"
      J. Venna, S. Kaski
    * "Learning a Parametric Embedding by Preserving Local Structure"
      L.J.P. van der Maaten

    Parameters
    ----------
    X : array, shape (n_samples, n_features) or (n_samples, n_samples)
        If the metric is 'precomputed' X must be a square distance
        matrix. Otherwise it contains a sample per row.

    X_embedded : array, shape (n_samples, n_components)
        Embedding of the training data in low-dimensional space.

    n_neighbors : int, optional (default: 5)
        Number of neighbors k that will be considered.

    precomputed : bool, optional (default: False)
        Set this flag if X is a precomputed square distance matrix.

    Returns
    -------
    trustworthiness : float
        Trustworthiness of the low-dimensional embedding.
    """
    if precomputed:
        dist_X = X
    else:
        dist_X = pairwise_distances(X, squared=True)
    dist_X_embedded = pairwise_distances(X_embedded, squared=True)
    ind_X = np.argsort(dist_X, axis=1)
    ind_X_embedded = np.argsort(dist_X_embedded, axis=1)[:, 1:n_neighbors + 1]

    n_samples = X.shape[0]
    t = 0.0
    ranks = np.zeros(n_neighbors)
    for i in range(n_samples):
        for j in range(n_neighbors):
            ranks[j] = np.where(ind_X[i] == ind_X_embedded[i, j])[0][0]
        ranks -= n_neighbors
        t += np.sum(ranks[ranks > 0])
    t = 1.0 - t * (2.0 / (n_samples * n_neighbors *
                          (2.0 * n_samples - 3.0 * n_neighbors - 1.0)))
    return t

    def make_plot(self):
        #plot gets arguments
        dates, prices = self.cmod.arguments_plot(buyerField=self.argCH_plot())
        print(dates)
        print(prices)
        #creating plot
        dates = np.array(dates)#converting list
        prices = np.array(prices)#converting list
        fig, self.plotTK = plt.subplots()
        s = np.argsort(dates)#hang price to date
        f = np.argsort(prices)#hang price to date
        self.plotTK.plot_date(dates[s], prices[f], 'bo-')

        self.plotTK.xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
        self.plotTK.fmt_xdata = DateFormatter('%Y-%m-%d %H:%M:%S')
        fig.autofmt_xdate()
        #merge plot and tkinter
        self.canvas = FigureCanvasTkAgg(fig, self.cview.frames[view.AboutPage].leftFrame)
        self.canvas.show()
        self.canvas.get_tk_widget().pack(side=tk.BOTTOM, fill=tk.BOTH, expand=True)
        #creating toolbar
        toolbar = NavigationToolbar2TkAgg(self.canvas, self.cview.frames[view.AboutPage].leftFrame)
        toolbar.update()
        #packing plot
        self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=True)

 def evaluation(self, test_data, test_label):
     dinx = np.array(list(self.train_drugs))
     DS = self.dsMat[:, dinx]
     tinx = np.array(list(self.train_targets))
     TS = self.tsMat[:, tinx]
     scores = []
     if self.K2 > 0:
         for d, t in test_data:
             if d in self.train_drugs:
                 if t in self.train_targets:
                     val = np.sum(self.U[d, :]*self.V[t, :])
                 else:
                     jj = np.argsort(TS[t, :])[::-1][:self.K2]
                     val = np.sum(self.U[d, :]*np.dot(TS[t, jj], self.V[tinx[jj], :]))/np.sum(TS[t, jj])
             else:
                 if t in self.train_targets:
                     ii = np.argsort(DS[d, :])[::-1][:self.K2]
                     val = np.sum(np.dot(DS[d, ii], self.U[dinx[ii], :])*self.V[t, :])/np.sum(DS[d, ii])
                 else:
                     ii = np.argsort(DS[d, :])[::-1][:self.K2]
                     jj = np.argsort(TS[t, :])[::-1][:self.K2]
                     v1 = DS[d, ii].dot(self.U[dinx[ii], :])/np.sum(DS[d, ii])
                     v2 = TS[t, jj].dot(self.V[tinx[jj], :])/np.sum(TS[t, jj])
                     val = np.sum(v1*v2)
             scores.append(np.exp(val)/(1+np.exp(val)))
     elif self.K2 == 0:
         for d, t in test_data:
             val = np.sum(self.U[d, :]*self.V[t, :])
             scores.append(np.exp(val)/(1+np.exp(val)))
     prec, rec, thr = precision_recall_curve(test_label, np.array(scores))
     aupr_val = auc(rec, prec)
     fpr, tpr, thr = roc_curve(test_label, np.array(scores))
     auc_val = auc(fpr, tpr)
     return aupr_val, auc_val

def scale_score(x, kind="quicksort", kind2="quicksort"):
    y = x.copy()
    order = np.argsort(x.flat, kind=kind)
    # Black magic ;-) Probably the smartest thing I came up with today.
    order_order = np.argsort(order, kind=kind2)
    y.flat[:] = order_order.astype(y.dtype)
    return y

def rowwise_rank(array, mask=None):
    """
    Take a 2D array and return the 0-indexed sorted position of each element in
    the array for each row.

    Example
    -------
    In [5]: data
    Out[5]:
    array([[-0.141, -1.103, -1.0171,  0.7812,  0.07  ],
           [ 0.926,  0.235, -0.7698,  1.4552,  0.2061],
           [ 1.579,  0.929, -0.557 ,  0.7896, -1.6279],
           [-1.362, -2.411, -1.4604,  1.4468, -0.1885],
           [ 1.272,  1.199, -3.2312, -0.5511, -1.9794]])

    In [7]: argsort(argsort(data))
    Out[7]:
    array([[2, 0, 1, 4, 3],
           [3, 2, 0, 4, 1],
           [4, 3, 1, 2, 0],
           [2, 0, 1, 4, 3],
           [4, 3, 0, 2, 1]])
    """
    # note that unlike scipy.stats.rankdata, the output here is 0-indexed, not
    # 1-indexed.
    return argsort(argsort(array))

 def _get_sorted_theta(self):
     '''sorts the integral points by bond in descending order'''
     depsf_arr = np.array([])
     V_f_arr = np.array([])
     E_f_arr = np.array([])
     xi_arr = np.array([])
     stat_weights_arr = np.array([])
     nu_r_arr = np.array([])
     r_arr = np.array([])
     for reinf in self.cont_reinf_lst:
         n_int = len(np.hstack((np.array([]), reinf.depsf_arr)))
         depsf_arr = np.hstack((depsf_arr, reinf.depsf_arr))
         V_f_arr = np.hstack((V_f_arr, np.repeat(reinf.V_f, n_int)))
         E_f_arr = np.hstack((E_f_arr, np.repeat(reinf.E_f, n_int)))
         xi_arr = np.hstack((xi_arr, np.repeat(reinf.xi, n_int)))
         stat_weights_arr = np.hstack((stat_weights_arr, reinf.stat_weights))
         nu_r_arr = np.hstack((nu_r_arr, reinf.nu_r))
         r_arr = np.hstack((r_arr, reinf.r_arr))
     argsort = np.argsort(depsf_arr)[::-1]
     # sorting the masks for the evaluation of F
     idxs = np.array([])
     for i, reinf in enumerate(self.cont_reinf_lst):
         idxs = np.hstack((idxs, i * np.ones_like(reinf.depsf_arr)))
     masks = []
     for i, reinf in enumerate(self.cont_reinf_lst):
         masks.append((idxs == i)[argsort])
     max_depsf = [np.max(reinf.depsf_arr) for reinf in self.cont_reinf_lst]
     masks = [masks[i] for i in np.argsort(max_depsf)[::-1]]
     return depsf_arr[argsort], V_f_arr[argsort], E_f_arr[argsort], \
             xi_arr[argsort], stat_weights_arr[argsort], \
             nu_r_arr[argsort], masks, r_arr[argsort]

def rforests(trainx, trainy, test, n_estimators=100, k=5):
	trainy = np.ravel(trainy)

	forest = RandomForestClassifier(n_estimators)
	forest.fit(trainx, trainy)


	prob_train = forest.predict_proba(trainx)
	prob_test = forest.predict_proba(test)

	# Since the index is the number of the country that's been chosen
	# we can use these with argsort to get the maximum 5., we will have to do this
	# for the entire matrix though.
	sort_train = np.argsort(prob_train)[:,-k:]
	sort_test = np.argsort(prob_test)[:,-k:]

	# Now we need to transform these back to countries, but to map I need to
	# have a dataframe.
	col_names = []

	for i in range(k):
		name = "country_destination_" + str(i+1)
		col_names.append(name)

	pred_train = pd.DataFrame(sort_train, columns=col_names)
	pred_test = pd.DataFrame(sort_test, columns=col_names)

	for name in col_names:
		pred_train[name] = pred_train[name].map(dicts.country)
		pred_test[name] = pred_test[name].map(dicts.country)

	pred_train = np.fliplr(pred_train)
	pred_test = np.fliplr(pred_test)

	return forest, pred_train, pred_test

    def show_heatmap(self, order_by = None,
                     order_by_row = None, order_by_col = None):
        if order_by:
            title = 'Network ordered by node covariate\n"%s"' % order_by
            o = np.argsort(self.node_covariates[order_by][:])
        elif order_by_row:
            title = 'Network ordered by row covariate\n"%s"' % order_by_row
            o = np.argsort(self.row_covariates[order_by_row][:])
        elif order_by_col:
            title = 'Network ordered by column covariate\n"%s"' % order_by_col
            o = np.argsort(self.col_covariates[order_by_col][:])
        else:
            title, o = 'Unordered adjacency matrix', np.arange(self.N)

        f, (ax_im, ax_ord) = plt.subplots(2, sharex = True)
        f.set_figwidth(3)
        f.set_figheight(6)
        A = self.adjacency_matrix()
        ax_im.imshow(A[o][:,o]).set_cmap('binary')
        ax_im.set_ylim(0, self.N - 1)
        ax_im.set_xticks([])
        ax_im.set_yticks([])
        ax_im.set_title(title)
        #plt.setp([ax_im.get_xticklabels(), ax_im.get_yticklabels()],
        #         visible = False)
        if order_by:
            ax_ord.scatter(np.arange(self.N), self.node_covariates[order_by][o])
            ax_ord.set_xlim(0, self.N - 1)
            ax_ord.set_ylim(self.node_covariates[order_by][o[0]],
                            self.node_covariates[order_by][o[-1]])
        plt.show()

def regenerate_dim(x):
    """ assume x in ns since epoch from the current time """
    msg = None  # msg allows us to see which shot/diag was at fault
    diffs = np.diff(x)
    # bincount needs a positive input and needs an array with N elts where N is the largest number input
    small = (diffs > 0) & (diffs < 1000000)
    sorted_diffs = np.sort(diffs[np.where(small)[0]])
    counts = np.bincount(sorted_diffs)
    bigcounts, bigvals = myhist(diffs[np.where(~small)[0]])

    if pyfusion.VERBOSE>0:
        print('[[diff, count],....]')
        print('small:', [[argc, counts[argc]] for argc in np.argsort(counts)[::-1][0:5]])
        print('big or negative:', [[bigvals[argc], bigcounts[argc]] for argc in np.argsort(bigcounts)[::-1][0:10]])

    dtns = 1 + np.argmax(counts[1:])  # skip the first position - it is 0
    # wgt0 = np.where(sorted_diffs > 0)[0]  # we are in ns, so no worry about rounding
    histo = plt.hist if pyfusion.DBG() > 1 else np.histogram
    cnts, vals = histo(x, bins=200)[0:2]
    # ignore the two end bins - hopefully there will be very few there
    wmin = np.where(cnts[1:-1] < np.max(cnts[1:-1]))[0]
    if len(wmin)>0:
        print('**********\n*********** Gap in data > {p:.2f}%'.format(p=100*len(wmin)/float(len(cnts))))
    x01111 = np.ones(len(x))  # x01111 will be all 1s except for the first elt.
    x01111[0] = 0
    errcnt = np.sum(bigcounts) + np.sum(np.sort(counts)[::-1][1:])
    if errcnt>0 or (pyfusion.VERBOSE > 0): 
        msg = str('** repaired length of {l:,}, dtns={dtns:,}, {e} erroneous utcs'
              .format(l=len(x01111), dtns=dtns, e=errcnt))

    fixedx = np.cumsum(x01111)*dtns
    wbad = np.where((x - fixedx)>1e8)[0]
    fixedx[wbad] = np.nan
    debug_(pyfusion.DEBUG, 3, key="repair", msg="repair of W7-X scrambled Langmuir timebase") 
    return(fixedx, msg)

	def target_neurons(self,nConnectPerInput,network,strCorr,bAntiCorr=False):
		numInput = self.dicProperties["IODim"]
		numNodesReservoir = self.dicProperties["ReservoirDim"]
		matTargetNeurons = np.zeros((numInput,nConnectPerInput))
		if strCorr == "Betweenness":
			self.lstBetweenness = betweenness_list(network)[0].a #get edge betweenness array
			lstSortedNodes = np.argsort(self.lstBetweenness)
			if not bAntiCorr:
				lstSortedNodes = lstSortedNodes[::-1]
			for i in range(numInput):
				lstRandIdx = rand_int_trunc_exp(0,numNodesReservoir,0.2,nConnectPerInput) # characteristic exponential decay is a fifth of the reservoir's size
				matTargetNeurons[i,:] = lstSortedNodes[lstRandIdx]
		elif "degree" in strCorr:
			# get the degree type
			idxDash = strCorr.find("-")
			strDegType = strCorr[:idxDash].lower()
			lstDegrees = degree_list(network,strDegType)
			# sort the nodes by their importance
			lstSortedNodes = np.argsort(lstDegrees)
			if not bAntiCorr:
				lstSortedNodes = lstSortedNodes[::-1]
			for i in range(numInput):
				lstRandIdx = rand_int_trunc_exp(0,numNodesReservoir,0.2,nConnectPerInput) # characteristic exponential decay is a fifth of the reservoir's size
				matTargetNeurons[i,:] = lstSortedNodes[lstRandIdx]
		else:
			matTargetNeurons = np.random.randint(0,numNodesReservoir,(numInput,nConnectPerInput))
		return matTargetNeurons.astype(int)

def argsort(x, topn=None, reverse=False):
    """Get indices of the `topn` smallest elements in array `x`.

    Parameters
    ----------
    x : array_like
        Array to sort.
    topn : int, optional
        Number of indices of the smallest(greatest) elements to be returned if given,
        otherwise - indices of all elements will be returned in ascending(descending) order.
    reverse : bool, optional
        If True - return the `topn` greatest elements, in descending order.

    Returns
    -------
    numpy.ndarray
        Array of `topn` indices that.sort the array in the required order.

    """
    x = np.asarray(x)  # unify code path for when `x` is not a np array (list, tuple...)
    if topn is None:
        topn = x.size
    if topn <= 0:
        return []
    if reverse:
        x = -x
    if topn >= x.size or not hasattr(np, 'argpartition'):
        return np.argsort(x)[:topn]
    # np >= 1.8 has a fast partial argsort, use that!
    most_extreme = np.argpartition(x, topn)[:topn]
    return most_extreme.take(np.argsort(x.take(most_extreme)))  # resort topn into order

    def SNfunc(self,data,sig,significancefloor=0.5):
        D=data.ravel()
        S=sig.ravel()

        args=numpy.argsort(-D/S)
        D=numpy.take(D,args)
        S=numpy.take(S,args)
        Dsum=numpy.cumsum(D)
        Ssum=numpy.cumsum(S**2)**0.5
        SN=(Dsum/Ssum).max()

        #regional SN
        import scipy.ndimage as  ndimage
        data[data/sig<significancefloor]=0
        masks, multiplicity = ndimage.measurements.label(data)
        labels=numpy.arange(1, multiplicity+1)
        SNs=numpy.zeros(multiplicity+1)
        SNs[0]=SN
        for i in range(multiplicity):
            D=data[masks==i+1].ravel()
            S=sig[masks==i+1].ravel()
            args=numpy.argsort(-D/S)
            D=numpy.take(D,args)
            S=numpy.take(S,args)
            Dsum=numpy.cumsum(D)
            Ssum=numpy.cumsum(S**2)**0.5
            SNi=(Dsum/Ssum).max()
            SNs[i+1]=SNi
        SNs=-numpy.sort(-SNs)
        return SNs

def get_heatmap(data_mat, name_for_saving_files,  pp,stimulus_on_time, stimulus_off_time,delta_ff, f0_start, f0_end):
    
    #Plot heatmap for validation 
    A1 = np.reshape(data_mat, (np.size(data_mat,0)*np.size(data_mat,1), np.size(data_mat,2)))
    if delta_ff == 1:
        delta_ff_A1 = np.zeros(np.shape(A1))
        for ii in xrange(0,np.size(A1,0)):
            delta_ff_A1[ii,:] = (A1[ii,:]-np.mean(A1[ii,f0_start:f0_end]))/(np.std(A1[ii,f0_start:f0_end])+0.1)
        B = np.argsort(np.mean(delta_ff_A1, axis=1))  
        print np.max(delta_ff_A1)
    else:
        B = np.argsort(np.mean(A1, axis=1)) 
        print np.max(A1)

    with sns.axes_style("white"):
        C = A1[B,:][-2000:,:]

        fig2 = plt.imshow(C,aspect='auto', cmap='jet', vmin = np.min(C), vmax = np.max(C))
        
        plot_vertical_lines_onset(stimulus_on_time)
        plot_vertical_lines_offset(stimulus_off_time)
        plt.title(name_for_saving_files)
        plt.colorbar()
        fig2 = plt.gcf()
        pp.savefig(fig2)
        plt.close()

    def __call__(self, filt, mask=None):
        '''
        Provide the iterator over the levels.
        '''
        self._check_filter(filt, mask)
        # This cover method is only for one-dimensional filter functions.
        assert(self.dim==1)
        # The interval length measures indices, not filter values
        # in this case.
        self.interval_length = 1. / \
            ( self.intervals[0] - (self.intervals[0]-1)*self.fract_overlap )
        self.step_size = self.interval_length*(1-self.fract_overlap)

        if mask is None:
            self.n = len(self.filt)
            self.sortorder = np.argsort(np.ravel(self.filt))
        else:
            idx = np.flatnonzero(mask)
            self.n = len(idx)
            sortorder = np.argsort(np.ravel(self.filt[mask]))
            self.sortorder = idx[sortorder]

        assert len(self.sortorder)==self.n

        self.iter = range(self.intervals[0]).__iter__()
        return self

 def predict_scores(self, test_data, N):
     dinx = np.array(list(self.train_drugs))
     DS = self.dsMat[:, dinx]
    # print DS drug-drug sim with 0 diagonal entries
     tinx = np.array(list(self.train_targets))
     TS = self.tsMat[:, tinx]
    # print TS target-target sim with 0 diagonal entries
     scores = []
     for d, t in test_data:
         if d in self.train_drugs: 
             if t in self.train_targets:
                 val = np.sum(self.U[d, :]*self.V[t, :])
             else:
                 jj = np.argsort(TS[t, :])[::-1][:N]
                 val = np.sum(self.U[d, :]*np.dot(TS[t, jj], self.V[tinx[jj], :]))/np.sum(TS[t, jj])
         else:
             if t in self.train_targets:
                 ii = np.argsort(DS[d, :])[::-1][:N]
                 val = np.sum(np.dot(DS[d, ii], self.U[dinx[ii], :])*self.V[t, :])/np.sum(DS[d, ii])
             else:
                 ii = np.argsort(DS[d, :])[::-1][:N]
                 jj = np.argsort(TS[t, :])[::-1][:N]
                 v1 = DS[d, ii].dot(self.U[dinx[ii], :])/np.sum(DS[d, ii])
                 v2 = TS[t, jj].dot(self.V[tinx[jj], :])/np.sum(TS[t, jj])
                 val = np.sum(v1*v2)
         if np.isnan(val):
             scores.append(0)
         else:
             scores.append(np.exp(val)/(1+np.exp(val)))
    # print smat #whole prediction matrix
     return np.array(scores) #from original code

def plotres(psr,deleted=False,group=None,**kwargs):
    """Plot residuals, compute unweighted rms residual."""

    res, t, errs = psr.residuals(), psr.toas(), psr.toaerrs
    
    if (not deleted) and N.any(psr.deleted != 0):
        res, t, errs = res[psr.deleted == 0], t[psr.deleted == 0], errs[psr.deleted == 0]
        print("Plotting {0}/{1} nondeleted points.".format(len(res),psr.nobs))

    meanres = math.sqrt(N.mean(res**2)) / 1e-6
    
    if group is None:
        i = N.argsort(t)
        P.errorbar(t[i],res[i]/1e-6,yerr=errs[i],fmt='x',**kwargs)
    else:
        if (not deleted) and N.any(psr.deleted):
            flagmask = psr.flagvals(group)[~psr.deleted]
        else:
            flagmask = psr.flagvals(group)

        unique = list(set(flagmask))
            
        for flagval in unique:
            f = (flagmask == flagval)
            flagres, flagt, flagerrs = res[f], t[f], errs[f]
            i = N.argsort(flagt)
            P.errorbar(flagt[i],flagres[i]/1e-6,yerr=flagerrs[i],fmt='x',**kwargs)
        
        P.legend(unique,numpoints=1,bbox_to_anchor=(1.1,1.1))

    P.xlabel('MJD'); P.ylabel('res [us]')
    P.title("{0} - rms res = {1:.2f} us".format(psr.name,meanres))