Normalize the eigenvectors of a reversible Markov state model according
    to our preferred scheme.
976081
def _normalize_eigensystem(u, lv, rv):
    
    # first normalize the stationary distribution separately
    lv[:, 0] = lv[:, 0] / np.sum(lv[:, 0])

    for i in range(1, lv.shape[1]):
        # the remaining left eigenvectors to satisfy
        # <\phi_i, \phi_i>_{\mu^{-1}} = 1
        lv[:, i] = lv[:, i] / np.sqrt(np.dot(lv[:, i], lv[:, i] / lv[:, 0]))

    for i in range(rv.shape[1]):
        # the right eigenvectors to satisfy <\phi_i, \psi_j> = \delta_{ij}
        rv[:, i] = rv[:, i] / np.dot(lv[:, i], rv[:, i])

    return u, lv, rv
138937
def normalize(v):
    ''''''
    if 0 == np.linalg.norm(v):
        return v
    return v / np.linalg.norm(v)
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
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]
592952
def _regroup(self):
        
        # clean up old maps
        for j in range(self.nout):
            self.inv_map[j] = []

        # find smallest divergence between input component i
        # and output component j of the cluster mixture density
        for i in range(self.nin):
            self.min_kl[i] = np.inf
            j_min = None
            for j in range(self.nout):
                kl = kullback_leibler(self.f.components[i], self.g.components[j])
                if kl < self.min_kl[i]:
                    self.min_kl[i] = kl
                    j_min = j
            assert j_min is not None
            self.inv_map[j_min].append(i)
364736
private HMM<O> iterate(int[][] sequences) {
        HMM<O> hmm = new HMM<>(numStates, numSymbols);
        hmm.symbols = symbols;

        // gamma[n] = gamma array associated to observation sequence n
        double gamma[][][] = new double[sequences.length][][];

        // a[i][j] = aijNum[i][j] / aijDen[i]
        // aijDen[i] = expected number of transitions from state i
        // aijNum[i][j] = expected number of transitions from state i to j
        double aijNum[][] = new double[numStates][numStates];
        double aijDen[] = new double[numStates];

        for (int k = 0; k < sequences.length; k++) {
            if (sequences[k].length <= 2) {
                throw new IllegalArgumentException(String.format("Traning sequence %d is too short.", k));
            }

            int[] o = sequences[k];
            double[][] alpha = new double[o.length][numStates];
            double[][] beta = new double[o.length][numStates];
            double[] scaling = new double[o.length];
            forward(o, alpha, scaling);
            backward(o, beta, scaling);

            double xi[][][] = estimateXi(o, alpha, beta);
            double g[][] = gamma[k] = estimateGamma(xi);

            int n = o.length - 1;
            for (int i = 0; i < numStates; i++) {
                for (int t = 0; t < n; t++) {
                    aijDen[i] += g[t][i];

                    for (int j = 0; j < numStates; j++) {
                        aijNum[i][j] += xi[t][i][j];
                    }
                }
            }
        }

        for (int i = 0; i < numStates; i++) {
            if (aijDen[i] == 0.0) // State i is not reachable
            {
                System.arraycopy(a[i], 0, hmm.a[i], 0, numStates);
            } else {
                for (int j = 0; j < numStates; j++) {
                    hmm.a[i][j] = aijNum[i][j] / aijDen[i];
                }
            }
        }

        /*
         * initial state probability computation
         */
        for (int j = 0; j < sequences.length; j++) {
            for (int i = 0; i < numStates; i++) {
                hmm.pi[i] += gamma[j][0][i];
            }
        }

        for (int i = 0; i < numStates; i++) {
            hmm.pi[i] /= sequences.length;
        }

        /*
         * emission probability computation
         */
        for (int i = 0; i < numStates; i++) {
            double sum = 0.0;

            for (int j = 0; j < sequences.length; j++) {
                int[] o = sequences[j];
                for (int t = 0; t < o.length; t++) {
                    hmm.b[i][o[t]] += gamma[j][t][i];
                    sum += gamma[j][t][i];
                }
            }

            for (int j = 0; j < numSymbols; j++) {
                hmm.b[i][j] /= sum;
            }
        }

        return hmm;
    }
867655
private List<SequenceState<State, Observation, Path>> computeViterbiSequence(
            List<TimeStep<State, Observation, Path>> timeSteps, int originalGpxEntriesCount,
            QueryGraph queryGraph) {
        final HmmProbabilities probabilities
                = new HmmProbabilities(measurementErrorSigma, transitionProbabilityBeta);
        final ViterbiAlgorithm<State, Observation, Path> viterbi = new ViterbiAlgorithm<>();

        logger.debug("\n=============== Paths ===============");
        int timeStepCounter = 0;
        TimeStep<State, Observation, Path> prevTimeStep = null;
        int i = 1;
        for (TimeStep<State, Observation, Path> timeStep : timeSteps) {
            logger.debug("\nPaths to time step {}", i++);
            computeEmissionProbabilities(timeStep, probabilities);

            if (prevTimeStep == null) {
                viterbi.startWithInitialObservation(timeStep.observation, timeStep.candidates,
                        timeStep.emissionLogProbabilities);
            } else {
                computeTransitionProbabilities(prevTimeStep, timeStep, probabilities, queryGraph);
                viterbi.nextStep(timeStep.observation, timeStep.candidates,
                        timeStep.emissionLogProbabilities, timeStep.transitionLogProbabilities,
                        timeStep.roadPaths);
            }
            if (viterbi.isBroken()) {
                String likelyReasonStr = "";
                if (prevTimeStep != null) {
                    Observation prevGPXE = prevTimeStep.observation;
                    Observation gpxe = timeStep.observation;
                    double dist = distanceCalc.calcDist(prevGPXE.getPoint().lat, prevGPXE.getPoint().lon,
                            gpxe.getPoint().lat, gpxe.getPoint().lon);
                    if (dist > 2000) {
                        likelyReasonStr = "Too long distance to previous measurement? "
                                + Math.round(dist) + "m, ";
                    }
                }

                throw new IllegalArgumentException("Sequence is broken for submitted track at time step "
                        + timeStepCounter + " (" + originalGpxEntriesCount + " points). "
                        + likelyReasonStr + "observation:" + timeStep.observation + ", "
                        + timeStep.candidates.size() + " candidates: "
                        + getSnappedCandidates(timeStep.candidates)
                        + ". If a match is expected consider increasing max_visited_nodes.");
            }

            timeStepCounter++;
            prevTimeStep = timeStep;
        }

        return viterbi.computeMostLikelySequence();
    }
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)
838795
def propose(self):
        

        # Locally store size of matrix
        dims = self.stochastic.value.shape

        # Add normal deviate to value and symmetrize
        dev = rnormal(
            0,
            self.adaptive_scale_factor *
            self.proposal_sd,
            size=dims)
        symmetrize(dev)

        # Replace
        self.stochastic.value = dev + self.stochastic.value
585739
def Hk(self, k, m_pred, P_pred):  # returns state iteration matrix
        

        return self.H[:, :, int(self.index[self.H_time_var_index, k])]
1045021
def rskew_normal(mu, tau, alpha, size=()):
    
    size_ = size or (1,)
    len_ = np.prod(size_)
    return flib.rskewnorm(
        len_, mu, tau, alpha, np.random.normal(size=2 * len_)).reshape(size)
61098
def rsdl_rn(self, AX, Y):
        """"""

        # The primal residual normalisation term is
        # max( ||A x^(k)||_2, ||B y^(k)||_2 ) and B = -(I I I ...)^T.
        # The scaling by sqrt(Nb) of the l2 norm of Y accounts for the
        # block replication introduced by multiplication by B
        return max((np.linalg.norm(AX), np.sqrt(self.Nb) * np.linalg.norm(Y)))
662045
public static void normalizeF( DMatrixRMaj A ) {
        double val = normF(A);

        if( val == 0 )
            return;

        int size = A.getNumElements();

        for( int i = 0; i < size; i++) {
            A.div(i , val);
        }
    }
393705
def summarize(self):
        

        doc = '''Markov state model
------------------
Lag time         : {lag_time}
Reversible type  : {reversible_type}
Ergodic cutoff   : {ergodic_cutoff}
Prior counts     : {prior_counts}

Number of states : {n_states}
Number of nonzero entries in counts matrix : {counts_nz} ({percent_counts_nz}%)
Nonzero counts matrix entries:
    Min.   : {cnz_min:.1f}
    1st Qu.: {cnz_1st:.1f}
    Median : {cnz_med:.1f}
    Mean   : {cnz_mean:.1f}
    3rd Qu.: {cnz_3rd:.1f}
    Max.   : {cnz_max:.1f}

Total transition counts :
    {cnz_sum} counts
Total transition counts / lag_time:
    {cnz_sum_per_lag} units
Timescales:
    [{ts}]  units
'''
        counts_nz = np.count_nonzero(self.countsmat_)
        cnz = self.countsmat_[np.nonzero(self.countsmat_)]

        return doc.format(
            lag_time=self.lag_time,
            reversible_type=self.reversible_type,
            ergodic_cutoff=self.ergodic_cutoff,
            prior_counts=self.prior_counts,
            n_states=self.n_states_,
            counts_nz=counts_nz,
            percent_counts_nz=(100 * counts_nz / self.countsmat_.size),
            cnz_min=np.min(cnz),
            cnz_1st=np.percentile(cnz, 25),
            cnz_med=np.percentile(cnz, 50),
            cnz_mean=np.mean(cnz),
            cnz_3rd=np.percentile(cnz, 75),
            cnz_max=np.max(cnz),
            cnz_sum=np.sum(cnz),
            cnz_sum_per_lag=np.sum(cnz)/self.lag_time,
            ts=', '.join(['{:.2f}'.format(t) for t in self.timescales_]),
            )
1169566
def _backward(self, x_dot_parameters):
        
        I, J, _ = self.x.shape
        return backward(self._lattice, x_dot_parameters, I, J,
                        self.state_machine.n_states)
318336
def _normalize(self):
        ""
        l,m,n = self.powers
        self.norm = np.sqrt(pow(2,2*(l+m+n)+1.5)*
                            pow(self.exponent,l+m+n+1.5)/
                            fact2(2*l-1)/fact2(2*m-1)/
                            fact2(2*n-1)/pow(np.pi,1.5))
        return
286041
def stationary_distribution(P, C=None, mincount_connectivity=0):
    
    # can be replaced by msmtools.analysis.stationary_distribution in next msmtools release
    from msmtools.analysis.dense.stationary_vector import stationary_distribution as msmstatdist
    if C is None:
        if is_connected(P, strong=True):
            return msmstatdist(P)
        else:
            raise ValueError('Computing stationary distribution for disconnected matrix. Need count matrix.')

    # disconnected sets
    n = np.shape(C)[0]
    ctot = np.sum(C)
    pi = np.zeros(n)
    # treat each weakly connected set separately
    sets = connected_sets(C, mincount_connectivity=mincount_connectivity, strong=False)
    for s in sets:
        # compute weight
        w = np.sum(C[s, :]) / ctot
        pi[s] = w * msmstatdist(P[s, :][:, s])
    # reinforce normalization
    pi /= np.sum(pi)
    return pi