Set the transition matrix according to self-loop probabilities.

        Parameters
        ----------
        p_self : None, float in (0, 1), or np.ndarray [shape=(n_labels,)]
            Optional self-loop probability(ies), used for Viterbi decoding
542360
def set_transition(self, p_self):
        '''
        '''
        if p_self is None:
            self.transition = None
        else:
            self.transition = np.empty((len(self._classes), 2, 2))
            if np.isscalar(p_self):
                self.transition = transition_loop(2, p_self)
            elif len(p_self) != len(self._classes):
                raise ParameterError('Invalid p_self.shape={} for vocabulary size={}'.format(p_self.shape, len(self._classes)))
            else:
                for i in range(len(self._classes)):
                    self.transition[i] = transition_loop(2, p_self[i])
6986
def _reset_internal_state(self):
        
        self.rotation = np.array([[-1., 0., 0.], [0., 1., 0.], [0., 0., -1.]])
        self.pos = np.zeros(3)  # (x, y, z)
        self.last_pos = np.zeros(3)
        self.grasp = False
712865
def _precompute_probabilities(self):
        

        d_graph = self.d_graph
        first_travel_done = set()

        nodes_generator = self.graph.nodes() if self.quiet \
            else tqdm(self.graph.nodes(), desc='Computing transition probabilities')

        for source in nodes_generator:

            # Init probabilities dict for first travel
            if self.PROBABILITIES_KEY not in d_graph[source]:
                d_graph[source][self.PROBABILITIES_KEY] = dict()

            for current_node in self.graph.neighbors(source):

                # Init probabilities dict
                if self.PROBABILITIES_KEY not in d_graph[current_node]:
                    d_graph[current_node][self.PROBABILITIES_KEY] = dict()

                unnormalized_weights = list()
                first_travel_weights = list()
                d_neighbors = list()

                # Calculate unnormalized weights
                for destination in self.graph.neighbors(current_node):

                    p = self.sampling_strategy[current_node].get(self.P_KEY,
                                                                 self.p) if current_node in self.sampling_strategy else self.p
                    q = self.sampling_strategy[current_node].get(self.Q_KEY,
                                                                 self.q) if current_node in self.sampling_strategy else self.q

                    if destination == source:  # Backwards probability
                        ss_weight = self.graph[current_node][destination].get(self.weight_key, 1) * 1 / p
                    elif destination in self.graph[source]:  # If the neighbor is connected to the source
                        ss_weight = self.graph[current_node][destination].get(self.weight_key, 1)
                    else:
                        ss_weight = self.graph[current_node][destination].get(self.weight_key, 1) * 1 / q

                    # Assign the unnormalized sampling strategy weight, normalize during random walk
                    unnormalized_weights.append(ss_weight)
                    if current_node not in first_travel_done:
                        first_travel_weights.append(self.graph[current_node][destination].get(self.weight_key, 1))
                    d_neighbors.append(destination)

                # Normalize
                unnormalized_weights = np.array(unnormalized_weights)
                d_graph[current_node][self.PROBABILITIES_KEY][
                    source] = unnormalized_weights / unnormalized_weights.sum()

                if current_node not in first_travel_done:
                    unnormalized_weights = np.array(first_travel_weights)
                    d_graph[current_node][self.FIRST_TRAVEL_KEY] = unnormalized_weights / unnormalized_weights.sum()
                    first_travel_done.add(current_node)

                # Save neighbors
                d_graph[current_node][self.NEIGHBORS_KEY] = d_neighbors
839179
def meval(self, f, *args, **kw):
        
        # !! PART OF ORIGINAL DOCSTRING INCOMPATIBLE WITH CLASS INTERFACE !!
        # Example
        # -------
        # We set up multiple stationary models, one for a reference (ground)
        # state, and two for biased states, and group them in a
        # MultiStationaryModel.
        # >>> from pyemma.thermo import StationaryModel, MEMM
        # >>> m_1 = StationaryModel(f=[1.0, 0], label='biased 1')
        # >>> m_2 = StationaryModel(f=[2.0, 0], label='biased 2')
        # >>> m_mult = MEMM([m_1, m_2], [0, 0], label='unbiased')
        # Compute the stationary distribution for the two biased models
        # >>> m_mult.meval('stationary_distribution')
        # [array([ 0.73105858,  0.26894142]), array([ 0.88079708,  0.11920292])]
        # We set up multiple Markov state models for different temperatures
        # and group them in a MultiStationaryModel.
        # >>> import numpy as np
        # >>> from pyemma.msm import MSM
        # >>> from pyemma.thermo import MEMM
        # >>> b = 20  # transition barrier in kJ / mol
        # >>> temps = np.arange(300, 500, 25)  # temperatures 300 to 500 K
        # >>> p_trans = [np.exp(- b / kT) for kT in 0.00831*temps ]
        # >>> # build MSMs for different temperatures
        # >>> msms = [MSM(P=np.array([[1.0-p, p], [p, 1.0-p]])) for p in p_trans]
        # >>> # build Multi-MSM
        # >>> msm_mult = MEMM(pi=msms[0].stationary_distribution, label='300 K', models=msms)
        # Compute the timescales and see how they decay with temperature
        # Greetings to Arrhenius.
        # >>> np.hstack(msm_mult.meval('timescales'))
        # array([ 1523.83827932,   821.88040004,   484.06386176,   305.87880068,
        #          204.64109413,   143.49286817,   104.62539128,    78.83331598])
        # !! END OF INCOMPATIBLE PART !!
        return [_call_member(M, f, *args, **kw) for M in self.models]
1002225
def _set_random_state_from(self, other):
        
        try:
            # this works if randgen is an instance of random.Random()
            self.randgen.setstate(other.randgen.getstate())
        except AttributeError:
            # this works if randgen is an instance of numpy.random.RandomState()
            self.randgen.set_state(other.randgen.get_state())

        return self
294392
def preprocess(options):

    assert options.bfile!=None, 'Please specify a bfile.'

    
    if options.compute_cov:
       assert options.bfile!=None, 'Please specify a bfile.'
       assert options.cfile is not None, 'Specify covariance matrix basename'
       print('Computing covariance matrix')
       t0 = time.time()
       computeCovarianceMatrix(options.plink_path,options.bfile,options.cfile,options.sim_type)
       t1 = time.time()
       print(('... finished in %s seconds'%(t1-t0)))
       print('Computing eigenvalue decomposition')
       t0 = time.time()
       eighCovarianceMatrix(options.cfile)
       t1 = time.time()
       print(('... finished in %s seconds'%(t1-t0)))

    """ computing principal components """
    if options.compute_PCs>0:
       assert options.ffile is not None, 'Specify fix effects basename for saving PCs'
       t0 = time.time()
       computePCs(options.plink_path,options.compute_PCs,options.bfile,options.ffile)
       t1 = time.time()
       print(('... finished in %s seconds'%(t1-t0)))

    """ fitting the null model """
    if options.fit_null:
        if options.nfile is None:
            options.nfile = os.path.split(options.bfile)[-1]
            warnings.warn('nfile not specifed, set to %s'%options.nfile)
        print('Fitting null model')
        assert options.pfile is not None, 'phenotype file needs to be specified'
        # read pheno
        Y = readPhenoFile(options.pfile,idx=options.trait_idx)
        # read covariance
        if options.cfile is None:
            cov = {'eval':None,'evec':None}
            warnings.warn('cfile not specifed, a one variance compoenent model will be considered')
        else:
            cov = readCovarianceMatrixFile(options.cfile,readCov=False)
            assert Y.shape[0]==cov['eval'].shape[0],  'dimension mismatch'
        # read covariates
        F = None
        if options.ffile is not None:
            F = readCovariatesFile(options.ffile)
            assert Y.shape[0]==F.shape[0], 'dimensions mismatch'
        t0 = time.time()
        fit_null(Y,cov['eval'],cov['evec'],options.nfile, F)
        t1 = time.time()
        print(('.. finished in %s seconds'%(t1-t0)))

    """ precomputing the windows """
    if options.precompute_windows:
        if options.wfile==None:
            options.wfile = os.path.split(options.bfile)[-1] + '.%d'%options.window_size
            warnings.warn('wfile not specifed, set to %s'%options.wfile)
        print('Precomputing windows')
        t0 = time.time()
        pos = readBimFile(options.bfile)
        nWnds,nSnps=splitGeno(pos,size=options.window_size,out_file=options.wfile+'.wnd')
        print(('Number of variants:',pos.shape[0]))
        print(('Number of windows:',nWnds))
        print(('Minimum number of snps:',nSnps.min()))
        print(('Maximum number of snps:',nSnps.max()))
        t1 = time.time()
        print(('.. finished in %s seconds'%(t1-t0)))

    # plot distribution of nSnps
    if options.plot_windows:
        print('Plotting ditribution of number of SNPs')
        plot_file = options.wfile+'.wnd.pdf'
        plt = pl.subplot(1,1,1)
        pl.hist(nSnps,30)
        pl.xlabel('Number of SNPs')
        pl.ylabel('Number of windows')
        pl.savefig(plot_file)
572414
def p(i, sample_size, weights):
        

        # Determine the initial pick values
        weight_i = weights[i]
        weights_sum = sum(weights)

        # Build a list of weights that don't contain the weight `i`. This list will
        # be used to build the possible picks before weight `i`.
        other_weights = list(weights)
        del other_weights[i]

        # Calculate the probability
        probability_of_i = 0
        for picks in range(0, sample_size):

            # Build the list of possible permutations for this pick in the sample
            permutations = list(itertools.permutations(other_weights, picks))

            # Calculate the probability for this permutation
            permutation_probabilities = []
            for permutation in permutations:

                # Calculate the probability for each pick in the permutation
                pick_probabilities = []
                pick_weight_sum = weights_sum

                for pick in permutation:
                    pick_probabilities.append(pick / pick_weight_sum)

                    # Each time we pick we update the sum of the weight the next
                    # pick is from.
                    pick_weight_sum -= pick

                # Add the probability of picking i as the last pick
                pick_probabilities += [weight_i / pick_weight_sum]

                # Multiply all the probabilities for the permutation together
                permutation_probability = reduce(
                    lambda x, y: x * y, pick_probabilities
                    )
                permutation_probabilities.append(permutation_probability)

            # Add together all the probabilities for all permutations together
            probability_of_i += sum(permutation_probabilities)

        return probability_of_i
320708
def gen_stochastic_matrix(size, random_state=None):
    
    if not type(size) is tuple:
        size = (1, size)

    assert len(size) == 2

    n = random_state.uniform(size=(size[0], size[1] - 1))
    n = np.concatenate([np.zeros((size[0], 1)), n, np.ones((size[0], 1))], axis=1)

    A = np.diff(np.sort(n))
    return A.squeeze()
1191386
def predict(self, parameters, viterbi):
        
        x_dot_parameters = np.einsum('ijk,kl->ijl', self.x, parameters)

        if not viterbi:
            alpha = forward_predict(self._lattice, x_dot_parameters,
                                    self.state_machine.n_states)
        else:
            alpha = forward_max_predict(self._lattice, x_dot_parameters,
                                        self.state_machine.n_states)

        I, J, _ = self.x.shape

        class_Z = {}
        Z = -np.inf

        for state, predicted_class in self.states_to_classes.items():
            weight = alpha[I - 1, J - 1, state]
            class_Z[self.states_to_classes[state]] = weight
            Z = np.logaddexp(Z, weight)

        return {label: np.exp(class_z - Z) for label, class_z in class_Z.items()}
96804
def _randDense(states, actions, mask):
    
    # definition of transition matrix : square stochastic matrix
    P = _np.zeros((actions, states, states))
    # definition of reward matrix (values between -1 and +1)
    R = _np.zeros((actions, states, states))
    for action in range(actions):
        for state in range(states):
            # create our own random mask if there is no user supplied one
            if mask is None:
                m = _np.random.random(states)
                r = _np.random.random()
                m[m <= r] = 0
                m[m > r] = 1
            elif mask.shape == (actions, states, states):
                m = mask[action][state]  # mask[action, state, :]
            else:
                m = mask[state]
            # Make sure that there is atleast one transition in each state
            if m.sum() == 0:
                m[_np.random.randint(0, states)] = 1
            P[action][state] = m * _np.random.random(states)
            P[action][state] = P[action][state] / P[action][state].sum()
            R[action][state] = (m * (2 * _np.random.random(states) -
                                _np.ones(states, dtype=int)))
    return(P, R)
893273
def transition_matrix_samples(self):
         """
        res = np.empty((self.nsamples, self.nstates, self.nstates), dtype=config.dtype)
        for i in range(self.nsamples):
            res[i, :, :] = self._sampled_hmms[i].transition_matrix
        return res
1130986
def __pred_for_np2d(self, mat, num_iteration, predict_type):
        """"""
        if len(mat.shape) != 2:
            raise ValueError('Input numpy.ndarray or list must be 2 dimensional')

        def inner_predict(mat, num_iteration, predict_type, preds=None):
            if mat.dtype == np.float32 or mat.dtype == np.float64:
                data = np.array(mat.reshape(mat.size), dtype=mat.dtype, copy=False)
            else:
                """change non-float data to float data, need to copy"""
                data = np.array(mat.reshape(mat.size), dtype=np.float32)
            ptr_data, type_ptr_data, _ = c_float_array(data)
            n_preds = self.__get_num_preds(num_iteration, mat.shape[0], predict_type)
            if preds is None:
                preds = np.zeros(n_preds, dtype=np.float64)
            elif len(preds.shape) != 1 or len(preds) != n_preds:
                raise ValueError("Wrong length of pre-allocated predict array")
            out_num_preds = ctypes.c_int64(0)
            _safe_call(_LIB.LGBM_BoosterPredictForMat(
                self.handle,
                ptr_data,
                ctypes.c_int(type_ptr_data),
                ctypes.c_int(mat.shape[0]),
                ctypes.c_int(mat.shape[1]),
                ctypes.c_int(C_API_IS_ROW_MAJOR),
                ctypes.c_int(predict_type),
                ctypes.c_int(num_iteration),
                c_str(self.pred_parameter),
                ctypes.byref(out_num_preds),
                preds.ctypes.data_as(ctypes.POINTER(ctypes.c_double))))
            if n_preds != out_num_preds.value:
                raise ValueError("Wrong length for predict results")
            return preds, mat.shape[0]

        nrow = mat.shape[0]
        if nrow > MAX_INT32:
            sections = np.arange(start=MAX_INT32, stop=nrow, step=MAX_INT32)
            # __get_num_preds() cannot work with nrow > MAX_INT32, so calculate overall number of predictions piecemeal
            n_preds = [self.__get_num_preds(num_iteration, i, predict_type) for i in np.diff([0] + list(sections) + [nrow])]
            n_preds_sections = np.array([0] + n_preds, dtype=np.intp).cumsum()
            preds = np.zeros(sum(n_preds), dtype=np.float64)
            for chunk, (start_idx_pred, end_idx_pred) in zip_(np.array_split(mat, sections),
                                                              zip_(n_preds_sections, n_preds_sections[1:])):
                # avoid memory consumption by arrays concatenation operations
                inner_predict(chunk, num_iteration, predict_type, preds[start_idx_pred:end_idx_pred])
            return preds, nrow
        else:
            return inner_predict(mat, num_iteration, predict_type)
935420
def directInitialMatrix(self):   
        
        
        #First initialize state codes and the mapping with states. 
        self.setStateCodes()  

        #For each state, calculate the indices of reached states and rates using the transition function.
        results  = imap(self.transitionStates, self.mapping.values())
        
         #Simpler alternative that uses less memory. 
         #Would be competitive if the conversion from dok to csr is faster.  
#        D = dok_matrix((self.size,self.size),dtype=float)
#        for index,(col,rate) in enumerate(results):
#            D.update({(index,c): r for c,r in zip(col,rate)})
#        return D.tocsr()

        #preallocate memory for the rows, cols and rates of the sparse matrix      
        rows = np.empty(self.size,dtype=int)
        cols = np.empty(self.size,dtype=int)
        rates = np.empty(self.size,dtype=float)
            
        #now fill the arrays with the results, increasing their size if current memory is too small.
        right = 0
        for index,(col,rate) in enumerate(results): #more robust alternative: in izip(self.mapping.keys(),results)
            left = right
            right += len(col)
            if right >= len(cols):
                new_capacity = int(round(right * 1.5))  #increase the allocated memory if the vectors turn out to be too small.
                cols.resize(new_capacity)
                rates.resize(new_capacity)
                rows.resize(new_capacity)
            rows[left:right] = index #since states are sorted, the index indeed corresponds to the state.
            cols[left:right] = col
            rates[left:right] = rate   
           
        #Place all data in a coo_matrix and convert to a csr_matrix for quick computations.
        return coo_matrix((rates[:right],(rows[:right],cols[:right])),shape=(self.size,self.size)).tocsr()
954894
def _to_p(self, mode):
        

        if self.mode.endswith("A"):
            chans = self.channels[:-1]
            alpha = self.channels[-1]
            self._secondary_mode = self.mode[:-1]
        else:
            chans = self.channels
            alpha = None
            self._secondary_mode = self.mode

        palette = []
        selfmask = chans[0].mask
        for chn in chans[1:]:
            selfmask = np.ma.mask_or(selfmask, chn.mask)
        new_chn = np.ma.zeros(self.shape, dtype=int)
        color_nb = 0

        for i in range(self.height):
            for j in range(self.width):
                current_col = tuple([chn[i, j] for chn in chans])
                try:
                    next(idx
                         for idx in range(len(palette))
                         if palette[idx] == current_col)
                except StopIteration:
                    idx = color_nb
                    palette.append(current_col)
                    color_nb = color_nb + 1

                new_chn[i, j] = idx

        if self.fill_value is not None:
            if self.mode.endswith("A"):
                current_col = tuple(self.fill_value[:-1])
                fill_alpha = [self.fill_value[-1]]
            else:
                current_col = tuple(self.fill_value)
                fill_alpha = []
            try:
                next(idx
                     for idx in range(len(palette))
                     if palette[idx] == current_col)
            except StopIteration:
                idx = color_nb
                palette.append(current_col)
                color_nb = color_nb + 1

            self.fill_value = [idx] + fill_alpha

        new_chn.mask = selfmask
        self.palette = palette
        if alpha is None:
            self.channels = [new_chn]
        else:
            self.channels = [new_chn, alpha]
        self.mode = mode
202025
def _update_inernal(self, varp_list):
        """"""
        # The probability for the binary variable for the same latent dimension of any of the models is on.
        if self.group_spike:
            self._b_prob_all = 1.-param_to_array(varp_list[0].gamma_group)
            [np.multiply(self._b_prob_all, 1.-vp.gamma_group, self._b_prob_all) for vp in varp_list[1:]]
        else:
            self._b_prob_all = 1.-param_to_array(varp_list[0].binary_prob)
            [np.multiply(self._b_prob_all, 1.-vp.binary_prob, self._b_prob_all) for vp in varp_list[1:]]
559660
def p_pragma_assign(self, p):
        ''
        p[0] = Pragma(PragmaEntry(p[3], p[5], lineno=p.lineno(1)),
                      lineno=p.lineno(1))
        p.set_lineno(0, p.lineno(1))
1241877
def update(self, Pi, Tij):
         """
        from msmtools import analysis as msmana

        # update transition matrix by copy
        self._Tij = np.array(Tij)
        assert msmana.is_transition_matrix(self._Tij), 'Given transition matrix is not a stochastic matrix'
        assert self._Tij.shape[0] == self._nstates, 'Given transition matrix has unexpected number of states '
        # reset spectral decomposition
        self._spectral_decomp_available = False

        # check initial distribution
        assert np.all(Pi >= 0), 'Given initial distribution contains negative elements.'
        assert np.any(Pi > 0), 'Given initial distribution is zero'
        self._Pi = np.array(Pi) / np.sum(Pi)