def get_od_pair_index_not_in_dataset(self, O, D):
        """Return indices to O (D) from whose elements an od pair is not included in the travel data
        see unittest for an example
        """
        from numpy import unique1d, setdiff1d, zeros_like, logical_and, logical_or, where
        
        assert O.shape == D.shape
        
        id_attributes = self.get_id_attribute()
        max_id = max(O.max(), D.max(), id_attributes.max())
        digits = len(str(max_id)) + 1
        multiplier = 10 ** digits

        ODpair = O * multiplier + D
        idpair = id_attributes[:, 0] * multiplier + id_attributes[:, 1]
        missing_pairs = setdiff1d( unique1d(ODpair), unique1d(idpair) )

        results = zeros_like(D)
        for pair in missing_pairs:
            results += logical_and( O == pair // multiplier, D == pair % multiplier)
        
        results += logical_or(O < id_attributes[:, 0].min(), O > id_attributes[:, 0].max())
        results += logical_or(D < id_attributes[:, 1].min(), D > id_attributes[:, 1].max())
        
        return where(results)

def uniquerows(X, decimalplaces, return_frequencies=False):
    """ Returns array consisting of unique rows and a list of the number of 
    times each row appeared.  Uniqueness is subject to rounding to the 
    specified number of decimal places. """
    # TODO there could be a serious bug with hash collisions in this implementation
    # TODO update this function to use the integer packing as a hash function?
    if len(X.shape) == 1:
        return np.unique1d(X)
    hashvec = np.random.random(X.shape[1])
    rowhashes = np.dot(np.around(X, decimalplaces), hashvec)
    # only appears to be consistent up to 12 decimal places
    rowhashes = np.around(rowhashes, 12)

    #rowhashes.sort()
    uniqs, inds = np.unique1d(rowhashes, return_index=True)
    # return_index of unique1d does not always return the first one

    # TODO faster way of doing this??
    # maybe something like diff(where(diff(sort(rowhashes))))
    # slow but simple way iterating through all rows
    if return_frequencies:
        ol = dict()
        for i in rowhashes:
            if ol.has_key(i):
                ol[i] += 1 
            else:
                ol[i] = 1
        return X[inds, :], ol.values()
    else: 
        return X[inds, :]

def unique2d(arr,axis=0):
    """function to sort and eliminate replicate rows/columns of an array. Extension to numpy's unique()
   
   USAGE:
   arr  : 2d array at this stage
   axis : sort by axis = [0],1 ([rows],cols)
    """

    I=[]

    # check shape of arr
    rows,cols=arr.shape
    
    # to sort by cols, transpose and do the same as you would for rows
    if axis==1:
	arr=arr.T

    for k in range(0,cols):
	i,d=n.unique1d(arr[:,k],return_index=True)
	I=n.hstack((I,i))
    I=n.unique1d(I)
    I=I.tolist()
    arr_out = arr[I,:]

    if axis==1:
	arr_out=arr_out.T

    return arr_out

def get_ld_grid(photband, **kwargs):
    """
    Retrieve an interpolating grid for the LD coefficients
    
    Check outcome:
    
    #>>> bands = ['GENEVA.U', 'GENEVA.B', 'GENEVA.G', 'GENEVA.V']
    #>>> f_ld_grid = get_ld_grid(bands)
    #>>> ff = pyfits.open(_atmos['file'])
    #>>> all(ff['GENEVA.U'].data[257][2:]==f_ld_grid(ff['GENEVA.U'].data[257][0],ff['GENEVA.U'].data[257][1])[0:5])
    #True
    #>>> all(ff['GENEVA.G'].data[257][2:]==f_ld_grid(ff['GENEVA.G'].data[257][0],ff['GENEVA.G'].data[257][1])[10:15])
    #True
    #>>> ff.close()
    
    #Make some plots:
    
    #>>> photband = ['GENEVA.V']
    #>>> f_ld = get_ld_grid(photband)
    #>>> logg = 4.0
    #>>> mu = linspace(0,1,100)
    #>>> p = figure()
    #>>> p = gcf().canvas.set_window_title('test of function <get_ld_grid>')
    #>>> for teff in linspace(9000,12000,19):
    #...    out = f_ld(teff,logg)
    #...    a1x,a2x,a3x,a4x, I_x1 = out.reshape((len(photband),5)).T
    #...    p = subplot(221);p = title('Interpolation of absolute intensities')
    #...    p = plot(teff,I_x1,'ko')
    #...    p = subplot(222);p = title('Interpolation of LD coefficients')
    #...    p = scatter(4*[teff],[a1x,a2x,a3x,a4x],c=range(4),vmin=0,vmax=3,cmap=cm.spectral,edgecolors='none')
    #...    p = subplot(223);p = title('Without absolute intensity')
    #...    p = plot(mu,ld_eval(mu,[a1x,a2x,a3x,a4x]),'-')
    #...    p = subplot(224);p = title('With absolute intensity')
    #...    p = plot(mu,I_x1*ld_eval(mu,[a1x,a2x,a3x,a4x]),'-')    
    
    """
    # -- retrieve the grid points (unique values)
    teffs, loggs = get_ld_grid_dimensions(**kwargs)
    teffs_grid = np.sort(np.unique1d(teffs))
    loggs_grid = np.sort(np.unique1d(loggs))
    coeff_grid = np.zeros((len(teffs_grid), len(loggs_grid), 5 * len(photband)))

    # -- get the FITS-file containing the tables
    gridfile = get_file(**kwargs)
    # -- fill the grid
    ff = pyfits.open(gridfile)
    for pp, iband in enumerate(photband):
        teffs = ff[iband].data.field("Teff")
        loggs = ff[iband].data.field("logg")
        for ii, (iteff, ilogg) in enumerate(zip(teffs, loggs)):
            indext = np.searchsorted(teffs_grid, iteff)
            indexg = np.searchsorted(loggs_grid, ilogg)
            # -- array and list are added for backwards compatibility with some
            #   pyfits versions
            coeff_grid[indext, indexg, 5 * pp : 5 * (pp + 1)] = np.array(list(ff[iband].data[ii]))[2:]
    ff.close()
    # -- make an interpolating function
    f_ld_grid = InterpolatingFunction([teffs_grid, loggs_grid], coeff_grid)
    return f_ld_grid

 def unique1d(a,return_indices=False):
     """Replacement for numpy's unique1d"""
     import numpy
     if return_indices:
         indices,uniq = numpy.unique1d(a,True)
         return uniq,indices
     else:
         return numpy.unique1d(a)

def agroupby(*args, **kwds):
    """A groupby function which accepts and returns arrays.
    All passed arrays are expected to be one dimensional
    and of the same shape. All of the arrays are grouped by
    `key(arg[0])` and then returned.  The returned arrays will
    be two dimensional with each row corresponding to a group.
    The size of the first dimension is equal to the number of
    groups, and the size of the second dimension is equal the
    the size of the largest groups.  All smaller groups are
    padded with the value of the keyword argument `fill_value`."""
    keyfunc = kwds.get('key', lambda a: a)
    fill_val = kwds.get('fill_value', 0.0)
    args = [a.copy() for a in args]
    argsort = sorted(enumerate(args[0]), key=compose(keyfunc,itemgetter(1)))
    indexsort = [index for index, item in argsort]
    args = [a.take(indexsort) for a in args]
    # calculate groups
    g_mask = keyfunc(args[0])
    g_set = unique1d(g_mask)
    g_max = max([g_mask[g_mask==g].shape[0] for g in g_set])
    g_args = [fill_val * ones((len(g_set), g_max), dtype=a.dtype) for a in args]
    for gix, gval in enumerate(g_set):
        for ga, a in izip(g_args, args):
            b = a[g_mask==gval]
            ga[gix,:len(b)] = b
    return tuple(g_args)

    def __init__(self, dataset, predictand, maxSubsetSize=None):
        """Constructor for a tree from a dataset of regressors (that which we split) and a predictand (that which we try to purify in the leaves).

        :type dataset: titus.producer.cart.Dataset
        :param dataset: dataset of regressors only
        :type predictand: 1-d Numpy array
        :param predictand: predictands in a separate array with the same number of rows as the ``dataset``
        :type maxSubsetSize: positive integer or ``None``
        :param maxSubsetSize: maximum size of subset splits of categorical regressors (approximation for optimization in ``categoricalEntropyGainTerm`` and ``categoricalNVarianceGainTerm``)
        """

        self.dataset = dataset
        self.predictand = predictand
        self.maxSubsetSize = maxSubsetSize

        self.datasetSize = len(self.predictand.data)

        if self.predictand.tpe == numbers.Real:
            try:
                self.predictandUnique = numpy.unique(self.predictand.data)
            except TypeError:
                self.predictandUnique = numpy.unique1d(self.predictand.data)

        elif self.predictand.tpe == basestring:
            if self.datasetSize > 0:
                self.predictandDistribution = []
                for category in xrange(len(self.predictand.intToStr)):
                    frac = 1.0 * numpy.sum(self.predictand.data == category) / len(self.predictand.data)
                    self.predictandDistribution.append(frac)
            else:
                self.predictandDistribution = [0.0] * len(self.predictand.intToStr)

        else:
            raise RuntimeError

def main():
    experiments = [ '1kmf-sndr0h=50km', '1kmf-zs-no-05XP', '1kmf-zs-no-mm-05XP', '1kmf-zs-no-mm', '1kmf-z-no-snd', '1kmf-z-no-v2' ]

    grid = goshen_1km_grid(bounds=(slice(100, 180), slice(90, 170)))
    temp = goshen_1km_temporal(start=14400)

    obs_file_names = ['psu_straka_mesonet.pkl', 'ttu_sticknet.pkl', 'asos.pkl']
    all_obs = loadObs(obs_file_names, temp.getDatetimes(aslist=True), grid, grid.getWidthHeight())

    ens_obs = {}
    for exp in experiments:
        ens_obs[exp] = cPickle.load(open("cold_pool_obs_%s.pkl" % exp, 'r'))

    mm_ids = np.unique1d(all_obs['id'])

    for id in mm_ids:
        id_idxs = np.where(all_obs['id'] == id)
        for ob_var, ens_var in [('temp', 't'), ('dewp', 'td')]:

            pylab.figure()
            pylab.plot(all_obs['time'][id_idxs], 5 / 9. * (all_obs[ob_var][id_idxs] - 32), 'k-', label='Observed')

            for exp_name, exp in ens_obs.iteritems():
                pylab.plot(all_obs['time'][id_idxs], exp[ens_var][id_idxs] - 273.15, label=exp_name)

            pylab.xticks(temp.getEpochs(aslist=True), temp.getStrings("%H%M", aslist=True), rotation=30)
            pylab.xlim(temp.getEpochs(aslist=True)[0], temp.getEpochs(aslist=True)[-1])
            pylab.legend(loc=1)

            pylab.savefig("mm_timeseries_%s_%s.png" % (ens_var, id))
            pylab.close()        

    return

    def __init__(self, dataset, predictand, maxSubsetSize=None):
        """Constructor for a dataset of predictors only; the
        predictand is a distinct Dataset.Field, provided separately.

        maxSubsetSize is used to limit splitting of categorical
        predictors; see categoricalEntropyGainTerm and
        categoricalNVarianceGainTerm.
        """

        self.dataset = dataset
        self.predictand = predictand
        self.maxSubsetSize = maxSubsetSize

        self.datasetSize = len(self.predictand.data)

        if self.predictand.tpe == numbers.Real:
            try:
                self.predictandUnique = numpy.unique(self.predictand.data)
            except TypeError:
                self.predictandUnique = numpy.unique1d(self.predictand.data)

        elif self.predictand.tpe == basestring:
            if self.datasetSize > 0:
                self.predictandDistribution = []
                for category in xrange(len(self.predictand.intToStr)):
                    frac = 1.0 * numpy.sum(self.predictand.data == category) / len(self.predictand.data)
                    self.predictandDistribution.append(frac)
            else:
                self.predictandDistribution = [0.0] * len(self.predictand.intToStr)

        else:
            raise RuntimeError

def plotObservationsComposite(obs, map, title, file_name):
    pylab.clf()
    colors = [ 'r', 'g', 'b', 'c', 'm', '#660099', '#ff9900', '#006666' ]

    ob_ids = np.unique1d(obs['id'])
    ttu_label = False

    for ob_id in ob_ids:
        ob_idxs = np.where(obs['id'] == ob_id)[0]
        these_obs = obs[ob_idxs]
        ob_xs, ob_ys = map(these_obs['longitude'], these_obs['latitude'])

        if ob_id[0] == "P":
            ob_num = int(ob_id[1]) - 1
            pylab.plot(ob_xs, ob_ys, 'o', mfc=colors[ob_num], mec=colors[ob_num], ms=3, label="NSSL MM (%s)" % ob_id)
        else:
            if not ttu_label:
                label = "TTU Sticknet"
                ttu_label = True
            else:
                label = None

            pylab.plot(ob_xs[0], ob_ys[0], 'ko', ms=3, label=label)

    drawPolitical(map, scale_len=5)

    pylab.legend(loc=3, numpoints=1, prop={'size':'medium'})
    pylab.title(title)
    pylab.savefig(file_name)
    return

def labelmeanfilter_str(ys, x):
    # works also for string labels in ys, but requires 1D
    # from mailing list scipy-user 2009-02-11
    unil, unilinv = np.unique1d(ys, return_index=False, return_inverse=True)
    labelmeans = np.array(ndimage.mean(x, labels=unilinv, index=np.arange(np.max(unil)+1)))
    arr3 = labelmeans[unilinv]
    return arr3

def plot(times, rms_difference, legend_loc, y_label, y_lim, colors, styles, title, file_name):
#   exp_names = { "1km-control-mod-05XP":"MM + MWR05XP", "1km-control-no-mm":"No MM", "1km-control-mm":"MM" }
    exp_names = { "3km-control":r"5 dBZ, 3 m s$^{-1}$", "3km-control-adapt=0.80":r"RTPS $\alpha$ = 0.80", "3km-control-adapt=1.00":r"RTPS $\alpha$ = 1.00", 
        "3km-control-r0h=12km":r"$r_{0h}$ = 12 km", "3km-control-7dBZ,5ms":r'$\sigma_Z$ = 7 dBZ, $\sigma_{v_r}$ = 5 m s$^{-1}$' }
    pylab.figure()
    pylab.axes((0.1, 0.125, 0.85, 0.8))

    all_good = [] 

    for exp_name in sorted(rms_difference.keys()):
        good_idxs = np.where(~np.isnan(rms_difference[exp_name]))[0]

        name, radar= exp_name.split(':')
        if len(good_idxs) > 0:
            pylab.plot(times[good_idxs], rms_difference[exp_name][good_idxs], color=colors[exp_name], linestyle=styles[exp_name], label="%s (%s)" % (exp_names[name], radar))

        all_good.append(good_idxs)

    all_good_idxs = np.unique1d(np.concatenate(tuple(all_good)))

    pylab.plot([times.min(), times.max()], [0, 0], color='k', linestyle=':')

    pylab.xlabel(r"Time (UTC)", size='large')
    pylab.ylabel(y_label, size='large')
    pylab.xlim((times.min(), times.max()))
    pylab.ylim(y_lim)
    pylab.xticks(times[all_good_idxs], [ (datetime(2009, 6, 5, 18, 0, 0) + timedelta(seconds=int(t))).strftime("%H%M") for t in times[all_good_idxs] ], size='large', rotation=30)
    pylab.yticks(size='large')
    pylab.legend(loc=legend_loc, prop={'size':'small'})
    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

 def get_near(self, point, d=1):
     coords = np.array(point*self.__divisions/self.__sizes, int)
     return np.unique1d([
         i for co in itertools.product(*[
             range(max(0, c-d), min(div, c+d+1))
             for c, div in zip(coords, self.__divisions)
             ]) for i in self.__get_cell(co)
         ])

 def z_slab(self, bottom, top, allTimes=True):
     """remove all trajectories that are not in the slab
         defined by [bottom, top]"""
     if allTimes:
         selection = np.unique1d(np.where(
             np.bitwise_and(
                 self.positions[:,:,-1]>bottom,
                 self.positions[:,:,-1]<top
                 ))[1])
     else:
         selection = np.unique1d(np.where(
             np.bitwise_and(
                 self.positions[:,:,-1].max(axis=0)>bottom,
                 self.positions[:,:,-1].min(axis=0)<top
                 )))
     self.trajs = self.trajs[selection]
     self.positions = self.positions[:,selection]

def _stats(input, labels = None, index = None, do_sum2=False):
    '''returns count, sum, and optionally sum^2 by label'''

    def single_group(vals):
        if do_sum2:
            return vals.size, vals.sum(), (vals * vals.conjugate()).sum()
        else:
            return vals.size, vals.sum()
        
    if labels is None:
        return single_group(input)

    # ensure input and labels match sizes
    input, labels = numpy.broadcast_arrays(input, labels)

    if index is None:
        return single_group(input[labels > 0])

    if numpy.isscalar(index):
        return single_group(input[labels == index])

    # remap labels to unique integers if necessary, or if the largest
    # label is larger than the number of values.
    if ((not numpy.issubdtype(labels.dtype, numpy.int)) or 
        (labels.min() < 0) or (labels.max() > labels.size)):
        unique_labels, new_labels = numpy.unique1d(labels, return_inverse=True)

        counts = numpy.bincount(new_labels)
        sums = numpy.bincount(new_labels, weights=input.ravel())
        if do_sum2:
            sums2 = numpy.bincount(new_labels, weights=(input * input.conjugate()).ravel())

        idxs = numpy.searchsorted(unique_labels, index)
        # make all of idxs valid
        idxs[idxs >= unique_labels.size] = 0
        found = (unique_labels[idxs] == index)
    else:
        # labels are an integer type, and there aren't too many, so
        # call bincount directly.
        counts = numpy.bincount(labels.ravel())
        sums = numpy.bincount(labels.ravel(), weights=input.ravel())
        if do_sum2:
            sums2 = numpy.bincount(labels.ravel(), weights=(input * input.conjugate()).ravel())

        # make sure all index values are valid
        idxs = numpy.asanyarray(index, numpy.int).copy()
        found = (idxs >= 0) & (idxs < counts.size)
        idxs[~ found] = 0

    counts = counts[idxs]
    counts[~ found] = 0
    sums = sums[idxs]
    sums[~ found] = 0
    if not do_sum2:
        return (counts, sums)    
    sums2 = sums2[idxs]
    sums2[~ found] = 0
    return (counts, sums, sums2)

def plotSubplots(times, data, legend_loc, y_label, y_lim, colors, styles, exp_names, title, file_name):
    pylab.figure(figsize=(10, 8))
    pylab.subplots_adjust(left=0.075, right=0.95, top=0.925, bottom=0.1, wspace=0.175, hspace=0.275)

    radars = sorted(list(set([ name.split(":")[-1] for name in data.keys() ])))

    n_rows = 1
    n_cols = (len(radars) + 1) / 2

    lines = {}

    for idx, radar in enumerate(radars):
        all_good = []

        pylab.subplot(n_rows, n_cols, idx + 1)
        for exp in sorted([ e for e in data.keys() if e.split(":")[-1] == radar ]):
            good_idxs = np.where(~np.isnan(data[exp]))[0]
            if len(good_idxs) > 0:
                name = exp.split(":")[0]
                line = pylab.plot(times[good_idxs], data[exp][good_idxs], color=colors[exp], label=exp_names[name])
                lines[exp_names[name]] = line

            all_good.append(good_idxs)
        all_good_idxs = np.unique1d(np.concatenate(tuple(all_good)))
        pylab.plot([times.min(), times.max()], [0, 0], color='k', linestyle=':')
        pylab.axvline(14400, color='k', linestyle=':')

        pylab.xlabel(r"Time (UTC)", size='xx-large')
        pylab.ylabel(y_label, size='xx-large')
        pylab.xlim((times.min(), times.max()))
        pylab.ylim(y_lim)
        unique_times = np.sort(np.unique1d(times))
        pylab.xticks(unique_times[::2], [ (datetime(2009, 6, 5, 18, 0, 0) + timedelta(seconds=int(t))).strftime("%H%M") for t in unique_times ][::2], rotation=30, size='xx-large')
        pylab.yticks(size='xx-large')

#       pylab.title(radar)

    labels, line_objs = zip(*lines.items())
#   pylab.gcf().legend(line_objs, labels, 'lower right', prop={'size':'medium'})
    pylab.legend(line_objs, labels, loc=legend_loc, prop={'size':'medium'})

    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

def ismember_newer(totest, members):
    """
    A setmember1d, which works for totest arrays with duplicate values
    """
    uniques_in_test, rev_idx = np.unique1d(totest, return_inverse=True)
    uniques_in_members_mask = np.setmember1d(uniques_in_test, members)
    # Use this instead if members is not unique
    # uniques_in_members_mask = setmember1d(uniques_in_test, unique1d(members))
    return uniques_in_members_mask[rev_idx]

 def renumber_array(array):
     uniq, inds = np.unique1d(array, return_inverse=True)
     #res = np.zeros_like(res2.shape)
     newuniq = range(len(uniq))
     res = np.array([newuniq[i] for i in inds]).astype(np.int16)
     res.shape = array.shape
     #for i in range(len(uniq)):
         #print "Segment %s has %d pixels" % (i, np.where( res == i)[0].size)
     return res

def decimate(vals, decimation):
    B = numpy.unique1d(vals)
    N = WhereIs(vals, B)

    # n contains the count of the number of elements for each bin
    # B contains the index location for each element in the array vals
    # NOTE: This histogram function will become obsolete soon by the numpy people.  The behavior will change.
    (n, B) = numpy.histogram(vals, B, new=False)
    return(DataTruncation(n, decimation, len(B), N))

def volumetricsample(P, res):
    inds = occupiedlist.pointstovoxels(P, res)
    #np.random.shuffle(inds)
    # TODO accessing a private method should move this to occupiedlist, or make 
    # it part of the public api of occupiedlist?
    ids = occupiedlist._pack(inds)
    # This introduces a systematic error as inds_unique is just the first index encountered.
    ids_unique, inds_unique = np.unique1d(ids, return_index=True)
    return P[inds_unique]