r"""
    Determine the number of states from a set of discrete trajectories

    Parameters
    ----------
    dtrajs : list of int-arrays
        discrete trajectories
493462
def number_of_states(dtrajs):
    
    """
    # determine number of states n
    nmax = 0
    for dtraj in dtrajs:
        nmax = max(nmax, np.max(dtraj))
    # return number of states
    return nmax + 1
1053117
def generate_states(self, initial_condition=None, with_noise=True, stateseq=None):
        
        from pybasicbayes.util.stats import sample_discrete
        # Generate from the prior and raise exception if unstable
        T, K, n = self.T, self.num_states, self.D_latent
        A = self.trans_matrix

        # Initialize discrete state sequence
        dss = -1 * np.ones(T, dtype=np.int32) if stateseq is None else stateseq.astype(np.int32)
        assert dss.shape == (T,)
        gss = np.empty((T,n), dtype='double')

        if initial_condition is None:
            if dss[0] == -1:
                dss[0] = sample_discrete(self.pi_0)
            gss[0] = self.init_dynamics_distns[dss[0]].rvs()
        else:
            dss[0] = initial_condition[0]
            gss[0] = initial_condition[1]

        for t in range(1,T):
            # Sample discrete state given previous continuous state
            if with_noise:
                # Sample discre=te state from recurrent transition matrix
                if dss[t] == -1:
                    dss[t] = sample_discrete(A[dss[t-1], :])

                # Sample continuous state given current discrete state
                gss[t] = self.dynamics_distns[dss[t-1]].\
                    rvs(x=np.hstack((gss[t-1][None,:], self.inputs[t-1][None,:])),
                        return_xy=False)
            else:
                # Pick the most likely next discrete state and continuous state
                if dss[t] == -1:
                    dss[t] = np.argmax(A[dss[t-1], :])

                gss[t] = self.dynamics_distns[dss[t-1]]. \
                    predict(np.hstack((gss[t-1][None,:], self.inputs[t-1][None,:])))
            assert np.all(np.isfinite(gss[t])), "SLDS appears to be unstable!"

        self.stateseq = dss
        self.gaussian_states = gss
236702
def data(self) -> Sequence[Tuple[float, float]]:
        
        return [(angle, prob) for angle, prob in zip(self._rabi_angles,
                                                     self._excited_state_probs)]
501535
def IpaS(d, a, tt_instance=True):
    '''
    '''

    if d == 1:
        M = _np.array([[1, 0], [a, 1]]).reshape((1, 2, 2, 1), order='F')
    else:
        M = [None] * d
        M[0] = _np.zeros((1, 2, 2, 2))
        M[0][0, :, :, 0] = _np.array([[1, 0], [a, 1]])
        M[0][0, :, :, 1] = _np.array([[0, a], [0, 0]])
        for i in xrange(1, d - 1):
            M[i] = _np.zeros((2, 2, 2, 2))
            M[i][:, :, 0, 0] = _np.eye(2)
            M[i][:, :, 1, 0] = _np.array([[0, 0], [1, 0]])
            M[i][:, :, 1, 1] = _np.array([[0, 1], [0, 0]])
        M[d - 1] = _np.zeros((2, 2, 2, 1))
        M[d - 1][:, :, 0, 0] = _np.eye(2)
        M[d - 1][:, :, 1, 0] = _np.array([[0, 0], [1, 0]])
    if tt_instance:
        M = _matrix.matrix.from_list(M)
    return M
767366
def __d2f(disc_coords, v):
        
        return np.array([disc_coords[0] / v.shape[0],
                         disc_coords[1] / v.shape[1],
                         disc_coords[2] / v.shape[2]])
683683
def get_transitions(self, indexes, history_length=1):
        
        frame_batch_shape = (
                [indexes.shape[0]]
                + list(self.state_buffer.shape[1:-1])
                + [self.state_buffer.shape[-1] * history_length]
        )

        past_frame_buffer = np.zeros(frame_batch_shape, dtype=self.state_buffer.dtype)
        future_frame_buffer = np.zeros(frame_batch_shape, dtype=self.state_buffer.dtype)

        for buffer_idx, frame_idx in enumerate(indexes):
            past_frame_buffer[buffer_idx], future_frame_buffer[buffer_idx] = self.get_frame_with_future(
                frame_idx, history_length
            )

        actions = self.action_buffer[indexes]
        rewards = self.reward_buffer[indexes]
        dones = self.dones_buffer[indexes]

        data_dict = {
            'observations': past_frame_buffer,
            'actions': actions,
            'rewards': rewards,
            'observations_next': future_frame_buffer,
            'dones': dones,
        }

        for name in self.extra_data:
            data_dict[name] = self.extra_data[name][indexes]

        return data_dict
784494
def nonempty_set(C, mincount_connectivity=0):
    
    # truncate to states with at least one observed incoming or outgoing count.
    if mincount_connectivity > 0:
        C = C.copy()
        C[np.where(C < mincount_connectivity)] = 0
    return np.where(C.sum(axis=0) + C.sum(axis=1) > 0)[0]
668687
def _nanstd(array, axis=None, ddof=0):
    """"""

    if isinstance(axis, tuple):
        array = _move_tuple_axes_first(array, axis=axis)
        axis = 0
    return bottleneck.nanstd(array, axis=axis, ddof=ddof)
1245183
def num_time_steps(self):
    """"""

    num_time_steps = sum(t.num_time_steps for t in self.trajectories)
    return num_time_steps + self.num_completed_time_steps
789385
def visualize_qual_stats_dict_single(D, dest, title):
    

    # first find out how many cycles were there. it is going to be about 101 for
    # hiseq runs, and 251 in miseq runs, but these values may change from run to
    # run. although all lanes are expected to have the identical number of cycles
    # the following code makes sure that the number_of_cycles variable holds the
    # longest one if there is a variation between the number of cycles between
    # lanes
    number_of_cycles = 0
    for pair in ['1', '2']:
        if pair not in D:
            continue
        for tile in D[pair]:
            if len(D[pair][tile]['mean']) > number_of_cycles:
                number_of_cycles = len(D[pair][tile]['mean'])


    fig = plt.figure(figsize = (12, 8))

    plt.rcParams.update({'axes.linewidth' : 0.9})
    plt.rc('grid', color='0.50', linestyle='-', linewidth=0.1)

    all_tiles = {'1': {'mean': [0] * number_of_cycles, 'count': [0] * number_of_cycles},
                 '2': {'mean': [0] * number_of_cycles, 'count': [0] * number_of_cycles}
                }

    for i in range(0, number_of_cycles):
        means_p1 = []
        counts_p1 = []
        means_p2 = []
        counts_p2 = []
        for tile_id in D['1']:
            tile = D['1'][tile_id]
            means_p1.append(tile['mean'][i])
            counts_p1.append(tile['count'][i])
            if '2' in D and D['2']:
                tile = D['2'][tile_id]
                means_p2.append(tile['mean'][i])
                counts_p2.append(tile['count'][i])

        all_tiles['1']['mean'][i] = numpy.mean(means_p1)
        all_tiles['1']['count'][i] = sum(counts_p1)

        if '2' in D and D['2']:
            all_tiles['2']['mean'][i] = numpy.mean(means_p2)
            all_tiles['2']['count'][i] = sum(counts_p2)

    colors = cm.get_cmap('RdYlGn', lut=256)

    plt.grid(True)

    plt.subplots_adjust(left=0.02, bottom = 0.03, top = 0.95, right = 0.98)

    plt.xticks(list(range(number_of_cycles / 10, number_of_cycles, number_of_cycles / 10)), rotation=90, size='xx-small')
    plt.ylim(ymin = 0, ymax = 42)
    plt.xlim(xmin = 0, xmax = number_of_cycles - 1)
    plt.yticks(list(range(5, 41, 5)), size='xx-small')

    plt.fill_between(list(range(0, number_of_cycles)), [42 for _ in range(0, number_of_cycles)], y2 = 0, color = colors(0), alpha = 0.2)
    plt.plot(all_tiles['1']['mean'], color = 'orange', lw = 6)

    read_number_percent_dropdown = [42 * (x / all_tiles['1']['count'][0]) for x in all_tiles['1']['count']]
    if not len(set(read_number_percent_dropdown)) <= 1:
        plt.fill_between(list(range(0, number_of_cycles)), read_number_percent_dropdown, y2 = 0, color = 'black', alpha = 0.08)

        plt.text(5, 2.5, '%s' % big_number_pretty_print(all_tiles['1']['count'][0]), alpha=0.5)
    else:
        plt.text(5, 2.5, '%s' % big_number_pretty_print(all_tiles['1']['count'][0]), alpha=0.5)

    if '2' in all_tiles and all_tiles['2']:
        plt.plot(all_tiles['2']['mean'], color = 'purple', lw = 6)

    plt.figtext(0.5, 0.97, '%s' % (title), weight = 'black', size = 'xx-large', ha = 'center')

    try:
        plt.savefig(dest + '.tiff')
    except:
        plt.savefig(dest + '.png')

    return (all_tiles['1']['mean'], all_tiles['2']['mean'])
1210943
def init_huang(X, n_clusters, dissim, random_state):
    """"""
    n_attrs = X.shape[1]
    centroids = np.empty((n_clusters, n_attrs), dtype='object')
    # determine frequencies of attributes
    for iattr in range(n_attrs):
        freq = defaultdict(int)
        for curattr in X[:, iattr]:
            freq[curattr] += 1
        # Sample centroids using the probabilities of attributes.
        # (I assume that's what's meant in the Huang [1998] paper; it works,
        # at least)
        # Note: sampling using population in static list with as many choices
        # as frequency counts. Since the counts are small integers,
        # memory consumption is low.
        choices = [chc for chc, wght in freq.items() for _ in range(wght)]
        # So that we are consistent between Python versions,
        # each with different dict ordering.
        choices = sorted(choices)
        centroids[:, iattr] = random_state.choice(choices, n_clusters)
    # The previously chosen centroids could result in empty clusters,
    # so set centroid to closest point in X.
    for ik in range(n_clusters):
        ndx = np.argsort(dissim(X, centroids[ik]))
        # We want the centroid to be unique, if possible.
        while np.all(X[ndx[0]] == centroids, axis=1).any() and ndx.shape[0] > 1:
            ndx = np.delete(ndx, 0)
        centroids[ik] = X[ndx[0]]

    return centroids
211451
def gdaOnes(shape, dtype, numGhosts=1):
    
    res = GhostedDistArray(shape, dtype)
    res.setNumberOfGhosts(numGhosts)
    res[:] = 1
    return res
293746
def objectcount(data, key):
    """"""
    objkey = key.upper()
    return len(data.dt[objkey])
967897
def slamdunkGeneralStatsTable(self):
        


        headers = OrderedDict()
        headers['counted'] = {
            'title': '{} Counted'.format(config.read_count_prefix),
            'description': '# reads counted within 3\'UTRs ({})'.format(config.read_count_desc),
            'shared_key': 'read_count',
            'min': 0,
            'format': '{:,.2f}',
            'scale': 'YlGn',
            'modify': lambda x: float(x) * config.read_count_multiplier,
        }
        headers['retained'] = {
            'title': '{} Retained'.format(config.read_count_prefix),
            'description': '# retained reads after filtering ({})'.format(config.read_count_desc),
            'shared_key': 'read_count',
            'min': 0,
            'format': '{:,.2f}',
            'scale': 'YlGn',
            'modify': lambda x: float(x) * config.read_count_multiplier,
        }
        headers['mapped'] = {
            'title': '{} Mapped'.format(config.read_count_prefix),
            'description': '# mapped reads ({})'.format(config.read_count_desc),
            'shared_key': 'read_count',
            'min': 0,
            'format': '{:,.2f}',
            'scale': 'YlGn',
            'modify': lambda x: float(x) * config.read_count_multiplier,
        }
        headers['sequenced'] = {
            'title': '{} Sequenced'.format(config.read_count_prefix),
            'description': '# sequenced reads ({})'.format(config.read_count_desc),
            'shared_key': 'read_count',
            'min': 0,
            'format': '{:,.2f}',
            'scale': 'YlGn',
            'modify': lambda x: float(x) * config.read_count_multiplier,
        }

        self.general_stats_addcols(self.slamdunk_data, headers)
480883
def get_trajectories(self, indexes, rollout_length):
        
        # assert indexes.shape[0] > 1, "There must be multiple indexes supplied"
        assert rollout_length > 1, "Rollout length must be greater than 1"

        batch_indexes = (
                indexes.reshape(1, indexes.shape[0]) - np.arange(rollout_length - 1, -1, -1).reshape(rollout_length, 1)
        )

        return self.get_transitions(batch_indexes)
1167443
def batchsd(trace, batches=5):
    

    if len(np.shape(trace)) > 1:

        dims = np.shape(trace)
        # ttrace = np.transpose(np.reshape(trace, (dims[0], sum(dims[1:]))))
        ttrace = np.transpose([t.ravel() for t in trace])

        return np.reshape([batchsd(t, batches) for t in ttrace], dims[1:])

    else:
        if batches == 1:
            return np.std(trace) / np.sqrt(len(trace))

        try:
            batched_traces = np.resize(trace, (batches, int(len(trace) / batches)))
        except ValueError:
            # If batches do not divide evenly, trim excess samples
            resid = len(trace) % batches
            batched_traces = np.resize(trace[:-resid], 
                (batches, len(trace[:-resid]) / batches))

        means = np.mean(batched_traces, 1)

        return np.std(means) / np.sqrt(batches)
127092
def k_modes(X, n_clusters, max_iter, dissim, init, n_init, verbose, random_state, n_jobs):
    """"""
    random_state = check_random_state(random_state)
    if sparse.issparse(X):
        raise TypeError("k-modes does not support sparse data.")

    X = check_array(X, dtype=None)

    # Convert the categorical values in X to integers for speed.
    # Based on the unique values in X, we can make a mapping to achieve this.
    X, enc_map = encode_features(X)

    n_points, n_attrs = X.shape
    assert n_clusters <= n_points, "Cannot have more clusters ({}) " \
                                   "than data points ({}).".format(n_clusters, n_points)

    # Are there more n_clusters than unique rows? Then set the unique
    # rows as initial values and skip iteration.
    unique = get_unique_rows(X)
    n_unique = unique.shape[0]
    if n_unique <= n_clusters:
        max_iter = 0
        n_init = 1
        n_clusters = n_unique
        init = unique

    results = []
    seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
    if n_jobs == 1:
        for init_no in range(n_init):
            results.append(k_modes_single(X, n_clusters, n_points, n_attrs, max_iter,
                                          dissim, init, init_no, verbose, seeds[init_no]))
    else:
        results = Parallel(n_jobs=n_jobs, verbose=0)(
            delayed(k_modes_single)(X, n_clusters, n_points, n_attrs, max_iter,
                                    dissim, init, init_no, verbose, seed)
            for init_no, seed in enumerate(seeds))
    all_centroids, all_labels, all_costs, all_n_iters = zip(*results)

    best = np.argmin(all_costs)
    if n_init > 1 and verbose:
        print("Best run was number {}".format(best + 1))

    return all_centroids[best], enc_map, all_labels[best], \
        all_costs[best], all_n_iters[best]