                         Changes in version 1.2                         

  - Added randomMarkovChain() (previously an unexported helper with a
    different interface), dirichletChain() and identityChain(), the
    generators random, dirichlet_process and identity of PyDTMC.
    randomMarkovChain() takes the exact number of zero probabilities
    (zeros) and a mask of probabilities fixed in advance;
    dirichletChain() draws every row from a truncated Dirichlet process
    (stick-breaking construction), with an optional diagonal bias, and
    accepts any positive concentration instead of PyDTMC's whole numbers
    between 1 and n. Both take a seed that makes a draw reproducible
    without touching the random number stream of the session; the
    distribution of their draws was checked against PyDTMC's.
  - Added timeCorrelations() and timeRelaxations(), the stationary time
    correlations of the observables defined by one or two observed
    sequences, and the expected value of such an observable after t
    steps from an initial distribution (PyDTMC's time_correlations and
    time_relaxations). They agree with a direct numpy computation
    to 1e-14. PyDTMC 9.0.0 returns wrong values for chains with complex
    eigenvalues as soon as a lag exceeds the number of states, because
    it then uses an eigendecomposition with eigenvectors that are not
    biorthonormal; here the powers are applied by repeated
    multiplication and squaring. Unlike PyDTMC, timeRelaxations() does
    not require a unique stationary distribution.
  - toFile() and fromFile() support the XML format of PyDTMC, so chains
    can be exchanged with it in both directions (checked on 40 chains);
    the chain's name is stored as an attribute that PyDTMC ignores.
    Writing uses only base R; reading requires the new suggested package
    xml2.
  - summary() for markovchain objects gains the argument details. With
    the default FALSE the output is unchanged; with TRUE a block of
    further properties follows (size, rank, classes, irreducibility,
    period, regularity, absorbing, reversible, stochastically monotone,
    symmetric, entropy rate, SLEM, spectral gap, Kemeny's constant),
    also returned in the details element of the result. On 40 test
    chains they agree with PyDTMC's, except Kemeny's constant, which
    PyDTMC defines differently (from the fundamental matrix of an
    absorbing chain) and which was checked against a direct eigenvalue
    computation instead.
  - Added toDot() and toMermaid(), which export the transition graph as
    Graphviz DOT or Mermaid text, with escaped state names and optional
    digits, minProbability and layout direction.
  - hittingProbabilities() gains the arguments solver, tol and maxIter
    (completes #203, whose targets half shipped earlier). On each target
    the states whose probability is not already decided by graph
    reachability need a linear system \eqn{(I - Q) h = R}, which up to
    this version was always solved by a doubled Neumann series: every
    step squared a dense matrix, at \eqn{O(m^3)} per step and up to 200
    steps. The new default, solver = "direct", is a single LU
    factorisation; solver = "bicgstab" is an unpreconditioned BiCGSTAB
    iteration on the sparse system, two sparse matrix-vector products
    per step, as asked in that report; solver = "doubling" keeps the old
    method, and is also the automatic fallback of "direct" on a
    numerically singular system. On a gambler's ruin chain with 1200
    states and one target, one call went from 22.4 to 0.10 seconds,
    and 3000 states became feasible at about one second. The default is
    also more accurate, not just faster: on the defining equation of the
    hitting probabilities the residual of "direct" is at machine
    precision (2e-16) and does not grow with the size of the chain,
    against 9e-14 and growing for the doubled series, and all three
    solvers were checked against the closed-form gambler's ruin
    probabilities. tol and maxIter, which were hard-coded at 1e-13
    and 200, are now arguments; note that an unpreconditioned BiCGSTAB
    can stall slightly above a very tight tol on an ill-conditioned
    chain, in which case the warning reports the residual actually
    reached. That warning now prints the residual in scientific
    notation, instead of showing every residual below 1e-6 as 0.000000.
    BiCGSTAB restarts from the current residual on a Lanczos breakdown
    or a stagnation, and switches to the direct solver with a warning if
    the breakdown persists.
  - Memory-safety audit of the C++ code (checked under AddressSanitizer,
    UndefinedBehaviorSanitizer and valgrind). Inputs that used to crash
    R, read past the end of a vector or loop forever now give an
    ordinary error: ctmcFit() with a data list that lacks the times or
    whose states and times differ in length;
    generatorToTransitionMatrix() with a non-square generator; the
    multinomial confidence intervals with counts not shaped like the
    transition matrix; markovchainFit(method = "map") with a hyperparam
    matrix without row or column names, or containing NA/NaN; a
    transition matrix containing NA/NaN (the stochasticity check never
    returned); an ictmc whose Q or range disagree with its states; and a
    markovchainList whose chains have inconsistent slots in
    rmarkovchain(parallel = TRUE). Valid inputs behave as before.
  - Correctness fixes found by the same audit: ctmcFit(byrow = FALSE)
    now scales the columns (not the rows) of the fitted matrix and
    returns a ctmc with byrow = FALSE, so the generator columns sum to
    zero; the stochasticity check no longer accepts a matrix with
    negative entries whose rows sum to 1 (a markovchain built from one
    now fails validation, as it should); the complex-number comparison
    used internally compared x^2 - y^2 instead of the squared modulus
    x^2 + y^2; and impreciseProbabilityatT() can be interrupted when a
    small error asks for many iterations.
  - The documentation of noofVisitsDist() now says what the function
    computes (#139): for each state, the expected fraction of the first
    N steps spent in it, (1/N) * sum_{k=1}^N P^k[i, ] (multiplied by N,
    the expected numbers of visits). Earlier versions described it as
    the joint pdf of the numbers of visits, which the package does not
    compute. The values returned are unchanged.
  - sanitize, in createSequenceMatrix() and markovchainFit(), now also
    accepts "uniform" and "absorbing" besides TRUE/FALSE (#213). The
    states that possibleStates adds but that are never observed as a
    departure state used to get, with sanitize = TRUE, a row of ones,
    i.e. a uniform distribution over every state: the chain was assumed
    to move from an unobserved state to anywhere with equal probability,
    which is an assumption about the data rather than a consequence of
    it, and the documentation ("put 1 in all rows having rowSum equal to
    zero") did not make that consequence visible. "absorbing" instead
    puts the whole mass on the diagonal, so those states are absorbing:
    the matrix is stochastic, as with TRUE, but no transition that was
    never observed is invented. "uniform" is an explicit synonym of
    TRUE, and TRUE/FALSE keep their exact previous meaning, so no
    existing call changes behaviour. In markovchainFit() the new mode
    works with every method, unlike the absorbingStates argument, which
    stays MLE-only as documented; where both apply they agree, because
    they share the same code path. The documentation of sanitize now
    describes what each value does to the row probabilities.
  - predict() and conditionalDistribution() for markovchain objects now
    stop with an informative message when the starting state is not one
    of the states of the chain (previously: an obscure 'names' attribute
    error), and predict() also checks newdata and n.ahead. The
    documentation of predict() now states how the prediction is made
    (most probable transition from the last element of newdata, ties
    broken at random) (#145).
  - Fixed several functions that ignored the byrow slot and read a
    column-stochastic matrix (byrow = FALSE) as if it were
    row-stochastic (completes the fix of #148, which had so far
    addressed the structural analysis only): markovchainSequence() and
    rmarkovchain() (both the C++ and the R implementation) simulated a
    different chain, and committorAB(), expectedRewards(),
    expectedRewardsBeforeHittingA(), firstPassage(),
    firstPassageMultiple(), noofVisitsDist(),
    is.stochasticallyMonotone() and autoLump() returned results for the
    transposed chain. lump() stopped with an error on column-stochastic
    chains because the stationary distributions were averaged along the
    wrong dimension. Results for chains with byrow = TRUE are unchanged;
    meanFirstPassageTime() and meanNumVisits() keep returning their
    output in the orientation of the chain, as documented.
  - Fixed fromFile() on XML files: probabilities written with 17 digits
    (as toFile() does) could come back one unit in the last place away
    from the value written on platforms where long double is no wider
    than double, such as arm64 macOS, because R's own number parser
    (as.numeric()) is not exact there. The probabilities are now
    converted by an exact decimal-to-double routine, so the round trip
    is exact on every platform (this also makes the reading of files
    written by PyDTMC independent of the platform). Nothing changes in
    the API.
  - Fixed ExpectedTime() on a continuous time Markov chain with two
    states: the matrix of the states other than the target has a single
    row and column, which was reduced to a number and made the function
    stop with an obscure error (invalid 'nrow' value).
  - Fixed thread safety in the RcppParallel paths of
    rmarkovchain(parallel = TRUE) on a markovchainList and of
    markovchainFit(method = "bootstrap", parallel = TRUE). Both workers
    called arma::sample(), which touches R's global RNG state
    (.Random.seed) and is not reentrant: concurrent calls from several
    threads introduced a data race that could silently bias the
    generated transitions. Each worker now uses a per-sequence
    std::mt19937_64 seeded from a single R::unif_rand() draw taken on
    the main thread, so set.seed() still controls reproducibility and
    the race is gone. As a side effect of the per-sequence seeding, the
    parallel path is now also bit-identical across runs for the same
    seed, regardless of how TBB splits the work. The sequence generation
    on a markovchainList also stops ignoring each chain's byrow slot: a
    column-stochastic chain used to be simulated with its matrix
    transposed, i.e. as a different chain -- the parallel analogue of
    the fix to the single-threaded path above. Internal
    micro-optimisations in the same two workers (cube held by reference,
    O(1) state lookup, pre-allocated row buffer, output stored
    contiguously) reduce the overhead of the parallel split; no public
    API changed.
  - markovchainFit() gains the argument num.cores, used only when method
    = "bootstrap" and parallel = TRUE. rmarkovchain(parallel = TRUE),
    markovchainFit(method = "bootstrap", parallel = TRUE) and the
    R-level helper .markovchainSequenceParallel() now all pick the
    thread count from the same resolver: an explicit num.cores wins,
    otherwise getOption("RcppParallel.numThreads"), getOption("Ncpus"),
    RCPP_PARALLEL_NUM_THREADS and OMP_NUM_THREADS are honoured in that
    order, and the fallback is min(2L, parallel::detectCores()) as CRAN
    policy requires. Previously rmarkovchain(parallel = TRUE) and the R
    socket cluster both defaulted to detectCores() - 1 while the
    parallel bootstrap used whatever thread count had been set earlier
    in the session; the three paths are now consistent.
  - Added redistribute(), propagating an initial distribution through
    steps transitions and returning the whole trajectory (a (steps + 1)
    x n matrix including the initial time, or only the last row with
    lastOnly = TRUE), as PyDTMC's mc_redistribute does. The initial
    distribution may be omitted (uniform), a single state name (point
    mass) or a numeric vector, named or not; column-stochastic chains
    are handled.
  - Added topologicalEntropy(), the logarithm of the Perron root of the
    chain's support (adjacency) matrix, in any base (default 2,
    consistent with entropyRate()), normalizedEntropyRate(), the entropy
    rate divided by the topological entropy (always in [0, 1]; 0 when
    either vanishes) and relaxationTime(), the reciprocal of
    spectralGap() (Inf for periodic chains). PyDTMC's mixing_rate is the
    first element of impliedTimescales() and is not duplicated.
  - autoplot.markovchain() gained a type argument: "graph" (the default,
    unchanged), "eigenvalues" (eigenvalues in the complex plane with the
    unit circle) and "flow" (evolution of the distribution, via
    redistribute(), with the new steps and initial arguments) and
    "comparison" (two or more chains on the same states, matched by
    name, side by side: transition-matrix heatmaps on a common scale, or
    with what = "stationary" their stationary distributions; new
    arguments other and what).
  - autoplot.markovchain() no longer triggers a ggplot2 warning about
    the unknown parameter linewidth when probabilities are shown on the
    edges: the border of the labels is now removed with linewidth
    (ggplot2 >= 3.5.0) or label.size (earlier versions). The graph no
    longer prints ggplot2's message "Scale for linewidth is already
    present" either: the linewidth scale is now added once instead of
    once per layer.
  - Added higherOrderLogLik() (closes #201), returning the
    log-likelihood, deviance, AIC and BIC of the model fitted by
    fitHigherOrder() so that models of different orders can be compared.
    It is the log-likelihood of the fitted mixture model, not
    necessarily a maximum: with the default method = "lsq",
    fitHigherOrder() chooses the weights by least squares on the
    stationary distribution, not by maximum likelihood (see method =
    "mle" below). A model of order k can only be evaluated from
    observation k + 1 on, so start lets several orders be evaluated on
    exactly the same observations, which is required for the criteria to
    be comparable. fitHigherOrder() itself and its return value are
    unchanged.
  - Added assessIndependence() (closes #109), the Anderson-Goodman
    (1957) chi-squared (method = "Pearson") or likelihood-ratio (method
    = "G") test of the null hypothesis that consecutive states are
    independent against a first-order Markov chain, i.e. a test of
    independence on the table of one-step transition counts. Only states
    actually observed as a departure or arrival state count towards the
    degrees of freedom; a degenerate table (e.g. a constant sequence)
    gives zero degrees of freedom and an NA p-value instead of an error.
    It complements verifyMarkovProperty() and assessOrder(), which test
    the first-order property against higher orders.
  - Fixed fitHigherOrder() returning the starting weights for orders 2
    and above. Its objective function subtracted order times the
    stationary distribution instead of once, which made it practically
    constant in the weights, so the optimizer never left its starting
    point (equal weights 1/order); and even with the objective corrected
    its value, of the order of 1e-7 on real sequences, is below the
    default tolerance of Rsolnp::solnp(), which stopped immediately. The
    objective is now the squared distance sum((sum_i lambda_i Q_i X -
    X)^2) divided by its value at the starting point. The weights
    returned for order >= 2 change (and are now the minimizers of the
    documented objective, checked against a brute-force grid); the
    signature and the structure of the returned list are unchanged.
  - fitHigherOrder() gains the argument method = c("lsq", "mle"). The
    default "lsq" is the existing least squares fit, so calls without
    the argument are unaffected. "mle" estimates the weights by maximum
    likelihood, keeping the same empirical lag matrices and the same
    returned list (lambda, Q, X): for fixed lag matrices the
    log-likelihood is concave in the weights, so the maximum is global,
    and it is computed by the EM algorithm for mixture weights, which
    does not need Rsolnp. On the Alofi rainfall and preproglucacon
    sequences the maximum likelihood weights raise the log-likelihood of
    orders two and three by up to about 11 and 28 units over the least
    squares ones (checked against a brute-force search over the simplex
    in the tests). This is not the mixture transition distribution model
    of Raftery (1985), whose single matrix shared by all lags is
    estimated together with the weights.
  - Added fitMTD(), maximum likelihood estimation of the mixture
    transition distribution (MTD) model of Raftery (1985), in which the
    next state depends on the last order states through a mixture, with
    non-negative lag weights, of a single transition matrix shared by
    all lags (r(r-1) + order - 1 parameters instead of r^order (r-1)).
    The likelihood is maximized by the EM algorithm of Lebre and
    Bourguignon (2008), implemented in C++ on the distinct observed
    patterns, fitting the orders 1, ..., order in turn so that each
    starts also from the previous one (the log-likelihood therefore
    never decreases with the order), with optional random restarts
    (nstart) against local maxima; start sets the observations the
    likelihood conditions on. It returns the weights, the matrix as a
    markovchain object, the log-likelihood with AIC and BIC, and the
    matrix in the layout of fitHigherOrder(), so the fit can be passed
    to higherOrderLogLik(), which now counts r(r-1) + order - 1
    parameters for MTD fits (unchanged for fitHigherOrder() fits). On
    the data of Berchtold and Raftery (2002) it reproduces all four MTD
    rows of Tables 2 and 3 (log-likelihood and BIC) and the estimates
    printed for the Koeberg MTD(2) model; this is part of the unit tests
    and of the higher_order_markov_chains vignette, which describes the
    model. Negative lag weights, allowed by Raftery's original
    formulation, are not supported.
  - Added higherOrderPredict() and higherOrderSimulate(): next-state
    probabilities given the most recent states, and simulation of a
    sequence, for the models fitted by fitHigherOrder() and fitMTD(),
    with the same transition probabilities that enter
    higherOrderLogLik() (checked in the tests: the sum of the log
    predicted probabilities equals the log-likelihood).
  - higherOrderLogLik() now accepts numeric and factor sequences,
    matched to the states of the fit as character strings, as fitMTD()
    does; previously it required a character vector, so
    higherOrderLogLik(x, fitMTD(x)) failed for a numeric x.
  - The documentation of fitHigherOrder() and higherOrderLogLik() now
    states that the least squares criterion is nearly flat in the
    weights (each lag matrix maps the empirical distribution onto itself
    up to end effects), that the "mle" weights are conditional on the
    empirical lag matrices, which in a mixture with free matrices would
    not make the weights identifiable.
  - For fits of fitHigherOrder(), higherOrderLogLik() now counts (r
    - 1)(1 + k(r - 1)) parameters, the dimension of the set of
    transition laws a mixture of k lag matrices can represent, instead
    of the naive k r(r - 1) + (k - 1), which exceeded it by r(k - 1)
    because the weights of such a mixture are not identifiable. The
    count is unchanged for k = 1 and for fits of fitMTD() (r(r - 1) + k
    - 1); AIC and BIC of orders two and above become smaller.
    higherOrderLogLik() is new in this version, so no released behaviour
    changes.
  - Fixed predictHommc(), which simulated from the wrong probabilities:
    it read the transition matrices by row although
    fitHighOrderMultivarMC() stores them by column (P[to, from]), so the
    draws used probabilities that did not sum to one and were silently
    renormalised by sample(); it used the past of the sequence being
    predicted instead of the past of each sequence that drives it; and
    it updated the state of each sequence before drawing the next ones
    in the same step. It now draws every sequence from the model of
    Ching et al. (2008) given the common past, and reads the matrices by
    row only for a hommc object with byrow = TRUE. The signature and the
    returned matrix are unchanged.
  - Fixed assessOrder() returning a statistic of 0 with 0 degrees of
    freedom, without any error, when the sequence is a factor: states
    are now compared as character strings, as in verifyMarkovProperty().
    The counting is vectorised (about 25 times faster on long
    sequences); results for character sequences are unchanged (checked
    on 200 random sequences).
  - markovchainFit() gains the argument progress = FALSE (closes #232).
    With progress = TRUE the "bootstrap" method shows a text progress
    bar (utils::txtProgressBar, style 3) over the simulation of the
    nboot bootstrap sequences and the estimation of a matrix from each
    of them, and can be interrupted while it runs. With parallel = TRUE
    the simulation runs in threads that cannot report progress, so the
    bar only covers the estimation step. The other methods are fast and
    ignore the argument. Results are unchanged (same random numbers with
    and without the bar), and the new argument is the last one, so
    existing calls are unaffected.
  - Added selectOrder(), selecting the order of a fully parameterized
    Markov chain by BIC (default) or AIC (Tong, 1975; Katz, 1981):
    chains of order 0, ..., maxOrder are fitted by maximum likelihood on
    the same observations (from start, by default maxOrder + 1), with
    the number of parameters counted as r^k (r - 1) or, with parameters
    = "observed", only for the transition probabilities not estimated as
    zero (the convention of Berchtold and Raftery, 2002). The table also
    reports the likelihood-ratio test of each order against the previous
    one. A list of sequences pools the counts without contexts crossing
    from one sequence to the next. It reproduces the independence and
    Markov chain rows of Tables 2 and 3 of Berchtold and Raftery (2002),
    and its order-1 test equals assessIndependence(method = "G"). Until
    now the package could only test order one against order two
    (assessOrder(), verifyMarkovProperty()).
  - The vignettes gain a section on comparing models of different orders
    with higherOrderLogLik() (higher_order_markov_chains) and one on
    assessIndependence() (an_introduction_to_markovchain_package), using
    the Alofi Island rainfall and preproglucacon sequences analysed by
    Avery and Henderson (1999).
  - Added to inst/extdata the Koeberg wind-direction and
    epileptic-seizure series of Berchtold and Raftery (2002), kindly
    provided by the authors (with a README giving the coding and the
    source). tests/testthat/test-berchtold-raftery-2002.R uses them to
    reproduce, up to the rounding of the published values, the
    log-likelihood and BIC of independence and of Markov chains of
    order 1-3 (Tables 2 and 3), the log-likelihood and BIC of the
    published MTD(2) model through higherOrderLogLik(), and the
    likelihood-ratio statistic of assessIndependence() implied by those
    tables. The published figures condition on the first 14
    observations, i.e. start = 15. Both vignettes show the reproduction.
  - Added toDictionary()/fromDictionary() and toFile()/fromFile(),
    converting a markovchain object to and from a plain R list, or a
    JSON/YAML/CSV file, mirroring PyDTMC's serialization functions of
    the same names. Unlike PyDTMC's own flat (from_state, to_state) ->
    probability dictionary ($n^2$ entries with no state grouping), the
    representation here nests the transition probabilities by source
    state (transitionMatrix[[from]][[to]]), which is both more compact
    to read and round-trips through R's own list-of-lists idiom without
    any special tuple-key handling. toDictionary() always normalizes to
    a row-stochastic, byrow = TRUE representation regardless of the
    source chain's own storage convention, so the dictionary is
    unambiguous on its own; fromDictionary() also accepts a plain matrix
    for transitionMatrix (honoring an explicit byrow, exactly as
    new("markovchain", ...) does) for building a dictionary by hand.
    JSON and YAML preserve the chain's name and every transition
    probability at full double precision (17 significant digits); CSV,
    which has no natural place for the chain's name, stores only the
    transition matrix -- the same limitation PyDTMC's own CSV format
    has. jsonlite/yaml are only in Suggests, and reading or writing the
    corresponding format raises an informative error if the package is
    not installed.
  - Added aggregateStates(), reducing a chain's state space to k
    macro-states by the spectral Kullback-Leibler minimization of Deng,
    Mehta and Meyn (2011), with both the spectral-bottom-up (grow from
    one macro-state) and spectral-top-down (merge down from every
    micro-state) heuristics, plus an adaptive choice between them
    following PyDTMC's own rule of thumb. Unlike PyDTMC, k is optional:
    left unset, it is chosen automatically via the eigengap of the
    transition matrix's eigenvalues, and the returned partition traces
    every macro-state back to its original state names (PyDTMC only
    labels the reduced states generically) so it can be fed straight
    into lump()/is.lumpable(). Verified numerically against lump(),
    including a subtlety worth flagging: the divergence is zero only
    when a source state's outgoing probability splits within a
    destination macro-state the same way regardless of source
    (proportional to the stationary distribution restricted to it) -- a
    stronger condition than the strong lumpability is.lumpable() checks,
    which only requires the macro-to-macro totals to agree. The
    classical Land of Oz chain is strongly lumpable and
    aggregateStates() correctly reproduces lump()'s aggregated chain for
    it, but its divergence is strictly positive for exactly this reason;
    this is documented, not a bug. The top-down heuristic's
    pairwise-merge candidate search was reimplemented as a worklist
    rather than a line-by-line port of PyDTMC's own version, which was
    found, by inspecting its source, to silently drop intermediate
    macro-states from a merge candidate whenever the two states chosen
    to merge were not adjacent in the current ordering.
  - Added birthDeath(), gamblersRuin() and urnModel(), three textbook
    Markov chain constructors: a general birth-death chain from
    adjacent-pair transition probabilities (with p/q deliberately of
    length n-1, so there is no leftover boundary probability to silently
    redistribute), the classical gambler's ruin (absorbing at both
    ends), and the Ehrenfest urn diffusion model (reflecting at both
    ends, exactly Binomial(balls, 0.5)-stationary).
  - Added toBoundedChain(), replacing a chain's first and last
    transition row with an absorbing, reflecting, or semi-reflecting
    rule, for use on any existing chain -- e.g. to add reflecting
    boundaries to a plain birthDeath() chain, which does not have them
    by construction.
  - Added tauchen() and rouwenhorst(), discretizing a stationary AR(1)
    process into a finite Markov chain by two classical econometric
    methods. Both return the actual numeric grid of the discretized
    process alongside the chain (not just generic state labels), since
    that is normally what the discretization was for. Verified
    numerically: rouwenhorst() reproduces the AR(1)'s theoretical
    unconditional variance and lag-1 autocorrelation to numerical
    precision at every chain size, including for highly persistent (near
    unit-root) processes where tauchen()'s fixed grid width is
    measurably less accurate.
  - Added sensitivity(), the closed-form sensitivity of the stationary
    distribution of a finite irreducible DTMC to a perturbation of one
    state's transition row, S[l,j] = pi[k]*(Z[l,j]-pi[j]) in terms of
    the fundamental matrix Z already used by kemenyConstant(). Unlike a
    positional, entry-by-entry sensitivity, S is built so that d %*% S
    gives the correct first-order change in pi for any zero-sum
    perturbation direction d of row k (Schweitzer, 1968; Meyer, 1980;
    Cho & Meyer, 2001), and the common two-state case (increase one
    transition, decrease another to compensate) is simply the difference
    of two rows of S.
  - Added mergeWith(), returning the convex combination (1-gamma)*P1 +
    gamma*P2 of two chains defined on the same state space. States are
    matched by name, not by matrix position: other's transition matrix
    is realigned to object's state order (and converted from
    column-stochastic if needed) before combining, and chains on
    different state-name sets are rejected with a clear error rather
    than silently combined positionally.
  - Documented the new sensitivity(), mergeWith(), birthDeath(),
    toBoundedChain(), urnModel(), gamblersRuin(), tauchen() and
    rouwenhorst() in the introductory vignette, with numerical examples
    checked against their original published sources: the closed-form
    gambler's ruin probability and the Ehrenfest chain's
    Binomial(balls, 0.5) stationary distribution (Grinstead &
    Snell, 2006), and a full replication of Kopecky and Suen's (2010)
    finding that rouwenhorst() matches a highly persistent AR(1)'s
    theoretical variance to numerical precision at every chain size
    while tauchen()'s fixed-width grid does not.
  - Added toNthOrder(), a thin, discoverability-only wrapper around the
    existing ^ operator (toNthOrder(object, order) is exactly object ^
    order), provided under this name because it is what PyDTMC calls the
    same operation; the underlying ^ method already uses repeated
    squaring, so there is nothing new to compute.
  - Added populationGeneticsModel(), building the Moran or Wright-Fisher
    population-genetics chain (loss/fixation of a mutant allele under
    mutation and selection) on n + 1 states. Verified against the
    classical neutral-drift fixation probability i/n for both models.
    Deliberately uses mutant relative fitness 1+s in both models, unlike
    PyDTMC's Moran implementation, which uses 1-s there while its own
    Wright-Fisher implementation uses 1+s -- an inconsistency this
    package avoids by using the same, more standard convention
    throughout; the two packages agree exactly at s = 0 and on the
    mutation-rate terms.
  - Added closestReversible() (closes #253), returning the reversible
    chain closest to a given one for a fixed stationary distribution,
    together with the distance it minimizes. It is the additive
    reversibilization (P + P*)/2 of Fill (1991), which is the orthogonal
    projection of P onto the reversible chains in the pi-weighted
    Hilbert-Schmidt norm of l^2(pi) and needs no numerical optimization:
    the stochasticity and non-negativity constraints hold automatically.
    The documentation states explicitly what is not solved -- the
    plain-Frobenius version of the same problem, and the harder problem
    of letting the stationary distribution vary (Nielsen & Weber, 2015).
  - markovchainFit() now accepts method = "laplace" for a list of
    sequences (closes #165), pooling the transition counts over the
    sequences and smoothing them exactly as it does for a single
    sequence. method = "bootstrap" on a list is still refused --
    resampling whole sequences and resampling transitions within them
    are different procedures with different standard errors, and the
    package will not pick one silently -- but the old "method not
    available for a list" is replaced by a message that says which
    methods to use instead.
  - Fixed markovchainFit(..., byrow = FALSE) returning an invalid
    markovchain object from every sequence/list fitting path: mle,
    laplace and map transposed the transition matrix but left the byrow
    slot at its TRUE default, while bootstrap did the mirror image (slot
    set to FALSE, matrix left row-stochastic). Such objects failed
    validObject() and were silently misread by every method that
    consults the slot (steadyStates(), is.irreducible(), and so on), so
    fits made with byrow = FALSE could yield wrong results downstream.
    The accompanying standardError and confidence-interval matrices now
    follow the estimate's orientation as well. Matrix/data.frame input
    is unchanged: there byrow describes the layout of the observations,
    and the fitted chain stays row-stochastic.
  - Added is.reversible(), checking detailed balance with respect to the
    stationary distribution of a finite irreducible DTMC. Only
    irreducibility is required, not aperiodicity, since reversibility
    and periodicity are independent properties.
  - Added mixingTime(), the total-variation mixing time of a finite
    irreducible, aperiodic DTMC. Unlike slem()/spectralGap(), this
    genuinely requires aperiodicity: a periodic chain's distribution
    never converges, so its mixing time is not defined and calling it on
    one raises a clear error instead of looping or returning a
    misleading value.
  - Added lazyChain(), building the lazy chain alpha*I + (1-alpha)*P; a
    standard device (Levin & Peres, 2017) for removing periodicity, e.g.
    to make mixingTime() applicable.
  - Added subchain() (closes #254), restricting a chain to a subset of
    states either as the raw principal submatrix (method = "submatrix",
    e.g. the Q block used by fundamentalMatrix()) or as a properly
    renormalized Markov chain conditioned on staying inside the subset
    (method = "renormalize"); the two are explicitly distinguished
    rather than conflated.
  - hittingProbabilities() gains the optional argument targets (part of
    #203): a character vector of state names for which the hitting
    probabilities are computed. Each column of the matrix is computed
    independently, so the cost is proportional to the number of targets
    instead of the number of states: on a random sparse chain of 300
    states, three targets take about 0.4 seconds instead of 45. The
    columns returned (rows for byrow = FALSE) are exactly those of the
    full matrix, in the order given; the default targets = NULL returns
    the full matrix as before, so existing calls and results are
    unaffected. A choice of numerical solvers for very large chains,
    also requested in #203, is not part of this change.

                        Changes in version 1.1.2                        

  - Added slem() and spectralGap() for finite irreducible DTMCs. Only
    irreducibility is required (not aperiodicity): a periodic chain
    correctly returns SLEM = 1 and spectral gap 0 instead of being
    rejected, since it genuinely does not contract towards its
    stationary distribution.
  - Added impliedTimescales(), returning the relaxation timescale of
    every non-trivial eigenvalue of a finite irreducible DTMC, with
    Inf/0 handled explicitly at the |lambda| = 1/0 boundaries.

                 Changes in version 1.1.1 (2026-09-18)                  

  - Added entropyRate() for finite irreducible DTMCs, with configurable
    logarithm base and support for row- and column-stochastic storage.
  - Added kemenyConstant() for finite irreducible DTMCs, with support
    for both row- and column-stochastic storage and tests for the
    zero-diagonal hitting-time convention.
  - Hardened public R and native C++ entry points against invalid
    indices, negative sizes, integer overflow, and unbounded
    simulations.

                        Changes in version 1.1.0                        

  - Added exact and approximate lumpability tools.
  - Made hittingProbabilities() robust to arbitrarily small positive
    transitions by combining graph-based zero/one classification with a
    relative-residual Neumann iteration.
  - Added fundamentalMatrix() for finite absorbing chains, including
    support for column-stochastic storage.
  - Improved absorptionProbabilities() by solving (I - Q) B = R directly
    instead of explicitly forming (I - Q)^{-1}.

                       Changes in version 0.11.1                        

  - Added explicit structural-zero support to assessStationarity().
  - Statistical inference functions now return R-standard htest objects
    and use the standard print.htest() output, while retaining the
    legacy dof field for compatibility.
  - Standardized assessOrder() on the same htest convention, including
    parameter = c(df = ...) and the standard test output.
  - Expanded the statistical_analysis section with numerical test
    results and time-homogeneity examples.
  - Added published literature examples for non-rejection of the Markov
    property and for empirical-versus-theoretical transition matrices.
  - Corrected structural-zero regression tests and the
    assessStationarity() documentation example used by package checks.
  - Corrected verifyEmpiricalToTheoretical() to use row totals when
    constructing expected transition counts for its row-wise multinomial
    likelihood-ratio test. This fixes a historical use of column totals
    and changes the Kullback et al. benchmark from 6.551795 to 6.518384,
    without changing its inferential conclusion.
  - Added explicit absorbingStates support to markovchainFit() for MLE
    fits of sequence data with known terminal states.

                       Changes in version 0.11.0                        

  - Fix for numerical issues
  - Removed deprecated .Dim methods
  - Add check monotonicity
  - Fixed a division-by-zero vulnerability in the C++ backend of
    generatorToTransitionMatrix(). The function now safely handles
    absorbing states (where the diagonal element is exactly 0) without
    returning NaNs or crashing.
  - Optimized memory allocation in the internal C++ function clean_nas()
    (used by createSequenceMatrix() and other fitting functions). By
    pre-calculating the exact number of valid elements before array
    instantiation, the package now avoids continuous memory
    reallocations, resulting in faster parsing of large sequences
    containing NA values.

                 Changes in version 0.10.3 (2026-02-02)                 

Handled R_CHECK_PACKAGES_USED_IN_DEMO

                       Changes in version 0.10.2                        

                       Changes in version 0.10.1                        

                 Changes in version 0.9.0 (2022-07-01)                  

  - Fixed CI calculations
