MarkovCorrelation.java

  1. package org.opentrafficsim.road.gtu.generator;

  2. import java.util.ArrayList;
  3. import java.util.LinkedHashMap;
  4. import java.util.LinkedHashSet;
  5. import java.util.List;
  6. import java.util.Locale;
  7. import java.util.Map;
  8. import java.util.Set;

  9. import org.djutils.exceptions.Throw;

  10. import nl.tudelft.simulation.jstats.streams.StreamInterface;

  11. /**
  12.  * Markov Chain functionality using state auto-correlations. Rather than specifying a Transition Matrix, this matrix is
  13.  * calculated from a steady-state and from given auto-correlations. Auto-correlation increases the probability that state S
  14.  * returns after state S. Correlation between states can be captured with grouping several states under a super state. This
  15.  * creates a Transition Matrix within a cell of a Transition Matrix. To the super matrix, this is simply a single
  16.  * state-supplying element, applicable to all previous states that are concerned within the element.<br>
  17.  * <br>
  18.  * This class is oblivious to intensity data of the states, in the sense that it must be provided to draw the next state. This
  19.  * class actually only remembers state auto-correlations. Together with input intensities, a part of the Transition Matrix is
  20.  * calculated on the fly. This flexibility over a fixed Transition Matrix allows 1) dynamic intensities, and 2) a single object
  21.  * of this class to be used for multiple processes in which intensities differ, but correlations are equal. This class is
  22.  * therefore also fairly flexible in terms of which states are concerned. Only states with correlation need to be added. For any
  23.  * state included in the input, for which no correlation is defined, a correlation of 0 is assumed. Reversely, states that are
  24.  * included in the class, but that are not part of the input states, are ignored.
  25.  * <p>
  26.  * Copyright (c) 2013-2020 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
  27.  * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
  28.  * <p>
  29.  * @version $Revision$, $LastChangedDate$, by $Author$, initial version 18 dec. 2017 <br>
  30.  * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
  31.  * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
  32.  * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
  33.  * @param <S> state type
  34.  * @param <I> intensity type
  35.  */
  36. public class MarkovCorrelation<S, I extends Number>
  37. {

  38.     /** Leaf node for each state present in the Markov Chain, including in all sub-groups. */
  39.     private Map<S, FixedState<S, I>> leaves = new LinkedHashMap<>();

  40.     /** Transition Matrix for each super state, i.e. the matrix within which states are put that have the given super state. */
  41.     private Map<S, TransitionMatrix<S, I>> superMatrices = new LinkedHashMap<>();

  42.     /** Matrix in which each state is located. */
  43.     private Map<S, TransitionMatrix<S, I>> containingMatrices = new LinkedHashMap<>();

  44.     /** Root matrix. */
  45.     private TransitionMatrix<S, I> root = new TransitionMatrix<>(null, 0.0); // input not important, it's for sub-groups

  46.     /**
  47.      * Adds a state to the root group of the Markov Chain. The correlation increases the chance that the state will occur after
  48.      * itself. We have: <br>
  49.      *
  50.      * <pre>
  51.      * p_ii = ss_i + (1 - ss_i) * c_i {eq. 1}
  52.      * </pre>
  53.      *
  54.      * where,
  55.      *
  56.      * <pre>
  57.      *   p_ii: the probability state i returns after state i
  58.      *   ss_i: the steady-state (overall mean) probability of state i
  59.      *   c_i:  correlation of state i
  60.      * </pre>
  61.      *
  62.      * Effective correlations of states depend on correlations of other states as well, so the given correlation is not
  63.      * guaranteed to result. One can easily see this from a system with 2 states: A and B. Suppose B has correlation. In order
  64.      * to maintain the same overall steady-state (occurrence proportion of A and B), it follows that state A also must be seen
  65.      * more frequently to follow itself.
  66.      *
  67.      * <pre>
  68.      * not correlated:  A B B A A A B A A A B A A B B
  69.      * very correlated: A A A A A A A A A B B B B B B (all B's grouped together, hence all A's grouped together and correlated)
  70.      * </pre>
  71.      *
  72.      * The effective correlation <i>c_i</i> of any state <i>i</i> can be calculated by reversing equation {@code eq. 1} using
  73.      * for <i>p_ii</i> the effective value after all correlations are applied. The procedure to derive the various probabilities
  74.      * from state <i>i</i> to state <i>j</i> (<i>p_ij</i>) is explained below. The procedure is based on the transition matrix
  75.      * <i>T</i>, in which each value gives the probability that the state changes from <i>i</i> (the row) to state <i>j</i> (the
  76.      * column). Consequently, the values in each row must sum to 1, as each state will be followed by another state.<br>
  77.      * <br>
  78.      * It is important that the transition matrix <i>T</i> results in a steady-state as provided. In particular we have for
  79.      * steady state <i>S</i> that <b><i>S</i>*<i>T</i> &#61; <i>S</i></b> should hold. Suppose we have <i>S</i> &#61; [0.7, 0.2,
  80.      * 0.1] for states A, B and C. Without any correlation this would give the base transition matrix:
  81.      *
  82.      * <pre>
  83.      *     | p_11  p_12  p_13 |   | 0.70  0.20  0.10 |
  84.      * T = | p_21  p_22  p_23 | = | 0.70  0.20  0.10 |
  85.      *     | p_31  p_32  p_33 |   | 0.70  0.20  0.10 |
  86.      * </pre>
  87.      *
  88.      * Our steady-state results as for whatever the previous state was, the steady-state probabilities are applied. Now suppose
  89.      * that state C has a correlation of 0.4. This would give that <i>p_33</i> :&#61; <i>p_33</i> + (1 - <i>p_33</i>) * <i>c</i>
  90.      * &#61; 0.46. With this increased value, the probabilities of row 3 no longer add up to 1. Hence, <i>p_31</i> and
  91.      * <i>p_32</i> should be reduced. However, we require that the same steady-state <i>S</i> is maintained. This will remain
  92.      * the case for as long as <i>T</i> remains a <i>reversible</i> Markov Chain. This means that each state has as much input
  93.      * probability, as it has output probability. A matrix where, except for the values on the diagonal, all column values are
  94.      * equal, is reversible. So the base <i>T</i> without correlation is reversible, and we only need to maintain reversibility.
  95.      * A method to maintain reversibility is to <i>scale symmetric pairs</i>. Hence, if we reduce <i>p_32</i>, we should reduce
  96.      * <i>p_23</i> by the same <i>factor</i>. Forcing row 3 to sum to 1, and scaling <i>p_31</i>, <i>p_13</i>, <i>p_32</i> and
  97.      * <i>p_23</i> by the same factor 0.6 we obtain the third matrix below.
  98.      *
  99.      * <pre>
  100.      *      | 0.70  0.20  0.10 |    | 0.70  0.20  0.10 |    | 0.70  0.20  0.06 |    | 0.74  0.20  0.06 |
  101.      * T =&gt; | 0.70  0.20  0.10 | =&gt; | 0.70  0.20  0.10 | =&gt; | 0.70  0.20  0.06 | =&gt; | 0.70  0.24  0.06 |
  102.      *      | 0.70  0.20  0.10 |    | 0.70  0.20  0.46 |    | 0.42  0.12  0.46 |    | 0.42  0.12  0.46 |
  103.      * </pre>
  104.      *
  105.      * As we reduce <i>p_13</i> and <i>p_23</i>, we also reduce the probability sums of rows 1 and 2. These reductions can be
  106.      * compensated by increasing the values on the diagonals, as is done in the fourth matrix. Note that changing the diagonal
  107.      * values does not affect reversibility. For example, 0.7*0.74 + 0.2*0.70 + 0.1*0.42 &#61; 0.7 for the first column.<br>
  108.      * <br>
  109.      * Changing the diagonal values <i>p_11</i> and <i>p_22</i> as the result of correlation for state C, shows that correlation
  110.      * of one state automatically introduces correlation at other states, as should also intuitively occur from the A-B example.
  111.      * The procedure can be started with the base <i>T</i> from steady-state <i>S</i> and can be repeated for each state
  112.      * <i>i</i> with correlation:
  113.      * <ol>
  114.      * <li>Increase <i>p_ii</i> to <i>p_ii</i> + (1-<i>p_ii</i>) * <i>c_i</i>.</li>
  115.      * <li>Reduce <i>p_ij</i> for all <i>j</i> unequal to <i>i</i> such that row <i>i</i> sums to 1, use one factor <i>f</i> for
  116.      * all.</li>
  117.      * <li>Reduce <i>p_ji</i> for all <i>j</i> unequal to <i>i</i> to maintain reversibility. Use factor <i>f</i> again.</li>
  118.      * <li>Set all <i>p_jj</i> for all <i>j</i> unequal to <i>i</i> such that row <i>j</i> sums to 1.</li>
  119.      * </ol>
  120.      * Knowing that each value <i>p_ij</i> gets reduced for correlation of state <i>i</i> and <i>j</i>, and realizing that the
  121.      * reduction factors <i>f</i> equal (1 - <i>c_i</i>) and (1 - <i>c_j</i>) respectively, the effective correlation can be
  122.      * calculated by adding all reductions in a row to the diagonal value <i>p_ii</i> and using {@code eq. 1}.<br>
  123.      * <br>
  124.      * See also "Construction of Transition Matrices of Reversible Markov Chains" by Qian Jiang.
  125.      * @param state S; state
  126.      * @param correlation double; correlation
  127.      * @throws IllegalArgumentException if correlation is not within the range (-1 ... 1), or the state is already defined
  128.      * @throws NullPointerException if state is null
  129.      */
  130.     public synchronized void addState(final S state, final double correlation)
  131.     {
  132.         Throw.whenNull(state, "State may not be null.");
  133.         Throw.when(this.leaves.containsKey(state), IllegalArgumentException.class, "State %s already defined.", state);
  134.         Throw.when(correlation <= -1.0 || correlation >= 1.0, IllegalArgumentException.class,
  135.                 "Correlation at root level need to be in the range (-1 ... 1).");
  136.         FixedState<S, I> node = new FixedState<>(state, correlation);
  137.         this.root.addNode(state, node);
  138.         this.containingMatrices.put(state, this.root);
  139.         this.leaves.put(state, node);
  140.     }

  141.     /**
  142.      * Adds a state to the group of the Markov Chain indicated by a super state. If the super state is not yet placed in a
  143.      * sub-group, the sub-group is created. Grouping is useful to let a set of states correlate to any other of the states in
  144.      * the set. For example, after state A, both states A and B can occur with some correlation, while state C is not correlated
  145.      * to states A and B. The same correlation is applied when the previous state was B, as it is also part of the same group.
  146.      * <br>
  147.      * <br>
  148.      * To explain sub-groups, suppose we have the following 3-state matrix in which the super state <i>s_2</i> is located (this
  149.      * can be the root matrix, or any sub-matrix). In this matrix, state <i>s_2</i> is replaced by a matrix <i>S_2</i>.
  150.      *
  151.      * <pre>
  152.      *       s_1   s_2   s_3             s_1   S_2   s_3
  153.      * s_1 | p_11  p_12  p_13 |    s_1 | p_11  p_12  p_13 |
  154.      * s_2 | p_21  p_22  p_23 | =&gt; S_2 | p_21  p_22  p_23 |
  155.      * s_3 | p_31  p_32  p_33 |    s_3 | p_31  p_32  p_33 |
  156.      * </pre>
  157.      *
  158.      * From the level of this matrix, nothing changes. Whenever the prior state was any of those inside <i>S_2</i>, row 2 is
  159.      * applied to determine the next state. If the next state is matrix <i>S_2</i>, the state is further specified by
  160.      * <i>S_2</i>. Matrix <i>S_2</i> itself will be:
  161.      *
  162.      * <pre>
  163.      *       s_2   s_4
  164.      * s_2 | p_22' p_24 |
  165.      * s_4 | p_42  p_44 |
  166.      * </pre>
  167.      *
  168.      * It will thus result in either state <i>s_2</i> or state <i>s_4</i>. More states can now be added to <i>S_2</i>, using the
  169.      * same super state <i>s_2</i>. In case the prior state was either <i>s_1</i> or <i>s_3</i>, i.e. no state included in the
  170.      * sub-group, the matrix of the sub-group defaults to fractions based on the steady-state only. Correlations are then also
  171.      * ignored.<br>
  172.      * <br>
  173.      * Correlation of <i>s_2</i> is applied to the whole group, and all other states of the group can be seen as sub-types of
  174.      * the group's super state. Correlations as considered inside the group (<i>c'_2</i> and <i>c'_4</i>) are mapped from the
  175.      * range <i>c_2</i> to 1. So, <i>c'_2</i> = 0, meaning that within the group there is no further correlation for the super
  176.      * state. Sub states, who are required to have an equal or larger correlation than the super state of the group, map
  177.      * linearly between <i>c_2</i> and 1.<br>
  178.      * <br>
  179.      * If the super state is only a virtual layer that should not in itself be a valid state of the system, it can simply be
  180.      * excluded from obtaining a new state using {@code getState()}.<br>
  181.      * <br>
  182.      * @param superState S; state of group
  183.      * @param state S; state to add
  184.      * @param correlation double; correlation
  185.      * @throws IllegalArgumentException if correlation is not within the range (0 ... 1), the state is already defined, or
  186.      *             superState is not yet a state
  187.      * @throws NullPointerException if an input is null
  188.      */
  189.     public synchronized void addState(final S superState, final S state, final double correlation)
  190.     {
  191.         Throw.whenNull(superState, "Super-state may not be null.");
  192.         Throw.whenNull(state, "State may not be null.");
  193.         Throw.when(this.leaves.containsKey(state), IllegalArgumentException.class, "State %s already defined.", state);
  194.         Throw.when(correlation < 0.0 || correlation >= 1.0, IllegalArgumentException.class,
  195.                 "Correlation at root level need to be in the range (-1 ... 1).");
  196.         if (!this.superMatrices.containsKey(superState))
  197.         {
  198.             // branch the super state in to a matrix
  199.             FixedState<S, I> superOriginal = this.leaves.get(superState);
  200.             // remove original from it's matrix
  201.             TransitionMatrix<S, I> superMatrix = this.containingMatrices.get(superState);
  202.             Throw.when(superMatrix == null, IllegalArgumentException.class, "No state has been defined for super-state %s.",
  203.                     superState);
  204.             superMatrix.removeNode(superState);
  205.             // replace with matrix
  206.             TransitionMatrix<S, I> matrix = new TransitionMatrix<>(superState, superOriginal.getCorrelation());
  207.             superMatrix.addNode(superState, matrix);
  208.             this.superMatrices.put(superState, matrix);
  209.             // add original node to that matrix
  210.             superOriginal.clearCorrelation();
  211.             matrix.addNode(superState, superOriginal);
  212.             this.containingMatrices.put(superState, matrix);
  213.         }
  214.         // add node
  215.         TransitionMatrix<S, I> superMatrix = this.superMatrices.get(superState);
  216.         Throw.when(correlation < superMatrix.getCorrelation(), IllegalArgumentException.class,
  217.                 "Sub states in a group can not have a lower correlation than the super state of the group.");
  218.         FixedState<S, I> node =
  219.                 new FixedState<>(state, (correlation - superMatrix.getCorrelation()) / (1.0 - superMatrix.getCorrelation()));
  220.         superMatrix.addNode(state, node);
  221.         this.containingMatrices.put(state, superMatrix);
  222.         this.leaves.put(state, node);
  223.         // register state as part of matrix node
  224.         this.root.registerInGroup(superState, state);
  225.         for (TransitionMatrix<S, I> matrix : this.superMatrices.values())
  226.         {
  227.             if (matrix.getState() != superState)
  228.             {
  229.                 matrix.registerInGroup(superState, state);
  230.             }
  231.         }
  232.     }

  233.     /**
  234.      * Draws a next state from this Markov Chain process, with predefined state correlations, but dynamic intensities. Any
  235.      * states that are present in the underlying Transition Matrix, but not present in the given states, are ignored. States
  236.      * that are not present in the underlying Transition Matrix, are added to it with a correlation of 0.
  237.      * @param previousState S; previous state
  238.      * @param states S[]; set of states to consider
  239.      * @param steadyState I[]; current steady-state intensities of the states
  240.      * @param stream StreamInterface; to draw random numbers
  241.      * @return S; next state
  242.      * @throws IllegalArgumentException if number of states is not the same as the stead-state length
  243.      * @throws NullPointerException if states, steadyState or stream is null
  244.      */
  245.     public synchronized S drawState(final S previousState, final S[] states, final I[] steadyState,
  246.             final StreamInterface stream)
  247.     {
  248.         Throw.whenNull(states, "States may not be null.");
  249.         Throw.whenNull(steadyState, "Steady-state may not be null.");
  250.         Throw.whenNull(stream, "Stream for random numbers may not be null.");
  251.         Throw.when(states.length != steadyState.length, IllegalArgumentException.class,
  252.                 "Number of states should match the length of the steady state.");
  253.         for (FixedState<S, I> node : this.leaves.values())
  254.         {
  255.             node.clearIntensity();
  256.         }
  257.         int n = states.length;
  258.         for (int i = 0; i < n; i++)
  259.         {
  260.             S state = states[i];
  261.             FixedState<S, I> leaf = this.leaves.get(state);
  262.             if (leaf == null)
  263.             {
  264.                 addState(state, 0.0);
  265.                 leaf = this.leaves.get(state);
  266.             }
  267.             leaf.setIntensity(steadyState[i]);
  268.         }
  269.         return this.root.drawState(previousState, stream);
  270.     }

  271.     /** {@inheritDoc} */
  272.     @Override
  273.     public String toString()
  274.     {
  275.         return "MarkovCorrelation [ " + this.root + " ]";
  276.     }

  277.     /**
  278.      * Base class for elements inside a Markov {@code TransitionMatrix}.
  279.      * <p>
  280.      * Copyright (c) 2013-2020 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
  281.      * <br>
  282.      * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
  283.      * <p>
  284.      * @version $Revision$, $LastChangedDate$, by $Author$, initial version 19 dec. 2017 <br>
  285.      * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
  286.      * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
  287.      * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
  288.      * @param <S> state type
  289.      * @param <I> intensity type
  290.      */
  291.     private abstract static class MarkovNode<S, I extends Number>
  292.     {

  293.         /** State. */
  294.         private final S state;

  295.         /** Correlation. */
  296.         private double correlation;

  297.         /**
  298.          * Constructor.
  299.          * @param state S; state
  300.          * @param correlation double; correlation
  301.          */
  302.         MarkovNode(final S state, final double correlation)
  303.         {
  304.             this.state = state;
  305.             this.correlation = correlation;
  306.         }

  307.         /**
  308.          * Returns the encapsulated state, which is either a fixed state, or the super-state of a group.
  309.          * @return S; encapsulated state, which is either a fixed state, or the super-state of a group
  310.          */
  311.         final S getState()
  312.         {
  313.             return this.state;
  314.         }

  315.         /**
  316.          * Returns the correlation.
  317.          * @return double; correlation
  318.          */
  319.         final double getCorrelation()
  320.         {
  321.             return this.correlation;
  322.         }

  323.         /**
  324.          * Clears the correlation.
  325.          */
  326.         protected final void clearCorrelation()
  327.         {
  328.             this.correlation = 0.0;
  329.         }

  330.         /**
  331.          * Returns the current intensity, used for the Markov Chain process.
  332.          * @return current intensity, used for the Markov Chain process, 0.0 if no intensity was provided
  333.          */
  334.         abstract double getIntensity();

  335.         /**
  336.          * Returns a state from this node, which is either a fixed state, or a randomly drawn state from a sub-group.
  337.          * @param previousState S; previous state
  338.          * @param stream StreamInterface; to draw random numbers
  339.          * @return S; state from this node, which is either a fixed state, or a randomly drawn state from a sub-group
  340.          */
  341.         abstract S drawState(S previousState, StreamInterface stream);

  342.     }

  343.     /**
  344.      * Transition matrix with functionality to return a next state, and to entail a set of fixed states mixed with matrices for
  345.      * sub-groups.
  346.      * <p>
  347.      * Copyright (c) 2013-2020 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
  348.      * <br>
  349.      * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
  350.      * <p>
  351.      * @version $Revision$, $LastChangedDate$, by $Author$, initial version 19 dec. 2017 <br>
  352.      * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
  353.      * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
  354.      * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
  355.      * @param <S> state type
  356.      * @param <I> intensity type
  357.      */
  358.     private static final class TransitionMatrix<S, I extends Number> extends MarkovNode<S, I>
  359.     {

  360.         /** List of state-defining states, where sub-groups are defined by a set of state. */
  361.         private List<Set<S>> states = new ArrayList<>();

  362.         /** List of nodes (fixed state or sub-group matrix). */
  363.         private List<MarkovNode<S, I>> nodes = new ArrayList<>();

  364.         /**
  365.          * Constructor.
  366.          * @param state S; super state representing the sub-group, or {@code null} for the root matrix.
  367.          * @param correlation double; correlation for the sub-group, or any value for the root matrix.
  368.          */
  369.         TransitionMatrix(final S state, final double correlation)
  370.         {
  371.             super(state, correlation);
  372.         }

  373.         /**
  374.          * Adds a node to the matrix.
  375.          * @param state S; state of the node
  376.          * @param node MarkovNode&lt;S, I&gt;; node
  377.          */
  378.         void addNode(final S state, final MarkovNode<S, I> node)
  379.         {
  380.             Set<S> set = new LinkedHashSet<>();
  381.             set.add(state);
  382.             this.states.add(set);
  383.             this.nodes.add(node);
  384.         }

  385.         /**
  386.          * Registers the state to belong to the group of super-state. This is used such that the matrix knows which previous
  387.          * states to map to the group.
  388.          * @param superState S; super-state of the group
  389.          * @param state S; state inside the group
  390.          */
  391.         void registerInGroup(final S superState, final S state)
  392.         {
  393.             for (Set<S> set : this.states)
  394.             {
  395.                 if (set.contains(superState))
  396.                 {
  397.                     set.add(state);
  398.                     return;
  399.                 }
  400.             }
  401.         }

  402.         /**
  403.          * Removes the node from the matrix, including group registration. This is used to replace a state with a group.
  404.          * @param state S; state to remove
  405.          */
  406.         void removeNode(final S state)
  407.         {
  408.             int i = -1;
  409.             for (int j = 0; j < this.states.size(); j++)
  410.             {
  411.                 if (this.states.get(j).contains(state))
  412.                 {
  413.                     i = j;
  414.                     break;
  415.                 }
  416.             }
  417.             if (i > -1)
  418.             {
  419.                 this.states.remove(i);
  420.                 this.nodes.remove(i);
  421.             }
  422.         }

  423.         /**
  424.          * Returns a state from this matrix. This is done by calculating the row of the Markov Transition Chain for the given
  425.          * previous state, and using those probabilities to draw an output state.
  426.          * @param previousState S; previous state
  427.          * @param stream StreamInterface; to draw random numbers
  428.          * @return S; state from this matrix
  429.          * @see MarkovCorrelation#addState(Object, double) algorithm for calculating the Transition Matrix
  430.          */
  431.         @Override
  432.         S drawState(final S previousState, final StreamInterface stream)
  433.         {
  434.             // figure out whether this matrix contains the previous state, and if so the correlation factor and row number i
  435.             boolean contains = false;
  436.             int n = this.states.size();
  437.             double intensitySum = 0.0;
  438.             int i = 0;
  439.             double iFactor = 1.0;
  440.             for (int j = 0; j < n; j++)
  441.             {
  442.                 intensitySum += this.nodes.get(j).getIntensity();
  443.                 if (this.states.get(j).contains(previousState))
  444.                 {
  445.                     i = j;
  446.                     contains = true;
  447.                     iFactor = 1.0 - this.nodes.get(i).getCorrelation();
  448.                 }
  449.             }
  450.             // gather probabilities and apply correlation factors
  451.             double[] p = new double[n];
  452.             double pSum = 0.0;
  453.             for (int j = 0; j < n; j++)
  454.             {
  455.                 if (i != j || !contains)
  456.                 {
  457.                     MarkovNode<S, I> node = this.nodes.get(j);
  458.                     double jFactor = 1.0;
  459.                     if (contains)
  460.                     {
  461.                         jFactor = 1.0 - node.getCorrelation();
  462.                     }
  463.                     p[j] = jFactor * iFactor * node.getIntensity() / intensitySum;
  464.                     pSum += p[j];
  465.                 }
  466.             }
  467.             // correct to get row sum = 1.0 by changing the diagonal value
  468.             if (contains)
  469.             {
  470.                 p[i] = 1.0 - pSum;
  471.             }
  472.             // make probabilities cumulative
  473.             for (int j = 1; j < n; j++)
  474.             {
  475.                 p[j] = p[j - 1] + p[j];
  476.             }
  477.             // draw
  478.             double r = stream.nextDouble();
  479.             for (int j = 0; j < n; j++)
  480.             {
  481.                 if (r < p[j])
  482.                 {
  483.                     return this.nodes.get(j).drawState(previousState, stream);
  484.                 }
  485.             }
  486.             throw new RuntimeException("Unexpected error while drawing state from matrix.");
  487.         }

  488.         /**
  489.          * Returns the current intensity, used for the Markov Chain process, as the sum of matrix elements.
  490.          * @return current intensity, used for the Markov Chain process, 0.0 if no intensity was provided
  491.          */
  492.         @Override
  493.         double getIntensity()
  494.         {
  495.             double intensity = 0;
  496.             for (MarkovNode<S, I> node : this.nodes)
  497.             {
  498.                 intensity += node.getIntensity();
  499.             }
  500.             return intensity;
  501.         }

  502.         /** {@inheritDoc} */
  503.         @Override
  504.         public String toString()
  505.         {
  506.             String superType = this.getState() == null ? "" : "(" + this.getState() + ")";
  507.             String statesStr = "";
  508.             String sep = "";
  509.             for (MarkovNode<S, I> node : this.nodes)
  510.             {
  511.                 statesStr += sep + node;
  512.                 sep = ", ";
  513.             }
  514.             return "T" + superType + "[ " + statesStr + " ]";
  515.         }
  516.     }

  517.     /**
  518.      * Container for a fixed state. Each state is reflected in a single object of this class. They are grouped in matrices,
  519.      * possibly all in the root. Subsets of all states may be grouped in a matrix, but no state is present in more than one
  520.      * matrix.
  521.      * <p>
  522.      * Copyright (c) 2013-2020 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
  523.      * <br>
  524.      * BSD-style license. See <a href="http://opentrafficsim.org/node/13">OpenTrafficSim License</a>.
  525.      * <p>
  526.      * @version $Revision$, $LastChangedDate$, by $Author$, initial version 19 dec. 2017 <br>
  527.      * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
  528.      * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
  529.      * @author <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
  530.      * @param <S> state type
  531.      * @param <I> intensity type
  532.      */
  533.     private static final class FixedState<S, I extends Number> extends MarkovNode<S, I>
  534.     {

  535.         /** Intensity. */
  536.         private I intensity;

  537.         /**
  538.          * Constructor.
  539.          * @param state S; state
  540.          * @param correlation double; correlation
  541.          */
  542.         FixedState(final S state, final double correlation)
  543.         {
  544.             super(state, correlation);
  545.         }

  546.         /**
  547.          * Returns the state.
  548.          * @param previousState S; previous state
  549.          * @param stream StreamInterface; to draw random numbers
  550.          * @return state S; the state
  551.          */
  552.         S drawState(final S previousState, final StreamInterface stream)
  553.         {
  554.             return getState();
  555.         }

  556.         /**
  557.          * Sets the current intensity.
  558.          * @param intensity I; intensity
  559.          */
  560.         void setIntensity(final I intensity)
  561.         {
  562.             this.intensity = intensity;
  563.         }

  564.         /**
  565.          * Clears the intensity.
  566.          */
  567.         void clearIntensity()
  568.         {
  569.             this.intensity = null;
  570.         }

  571.         /** {@inheritDoc} */
  572.         @Override
  573.         double getIntensity()
  574.         {
  575.             return this.intensity == null ? 0.0 : this.intensity.doubleValue();
  576.         }

  577.         /** {@inheritDoc} */
  578.         @Override
  579.         public String toString()
  580.         {
  581.             return String.format(Locale.US, "%s(%.2f)", getState(), getCorrelation());
  582.         }
  583.     }

  584. }