*/
public static void markovgen(final double[][] atrans, final int[] out, final int istart, final int seed) {
int i, ilo, ihi, ii, j, m = atrans.length, n = out.length;
double[][] cum = buildMatrix(atrans); // temporary matrix to hold cumulative probabilities
double r;
Ran ran = new Ran(seed); // use the random number generator Ran
if (m != atrans[0].length) throw new IllegalArgumentException("transition matrix must be square");
for (i=0; i<m; i++) { // fill cum and die if clearly not a transition matrix
for (j=1; j<m; j++) cum[i][j] += cum[i][j-1];
if (abs(cum[i][m-1]-1.) > 0.01)
throw new IllegalArgumentException("transition matrix rows must sum to 1");
}
j = istart; // the current state is kept in j
out[0] = j;
for (ii=1; ii<n; ii++) { // Μain loop
r = ran.doub()*cum[j][m-1]; // Slightly-off normalization gets corrected here
ilo = 0;
ihi = m;
while (ihi-ilo > 1) { // Use bisection to find location among the cumulative probabilities
i = (ihi+ilo) >> 1;
if (r>cum[j][i-1]) ilo = i;