From b01b0ab15bd254d27d68f08b839df310688bc751 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 15 Mar 2016 18:50:57 -0400 Subject: [PATCH 01/23] Add node labels to DiGraph --- quantecon/graph_tools.py | 74 +++++++++++++++++++++++++++++++++++----- 1 file changed, 66 insertions(+), 8 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index e75aeeecb..1fd68a9e0 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -27,6 +27,10 @@ class DiGraph(object): weighted : bool, optional(default=False) Whether to treat `adj_matrix` as a weighted adjacency matrix. + node_labels : array_like(ndim=1, default=None) + Array_like of length n containing the label associated with each + node. If None, the labels default to integers 0 through n-1. + Attributes ---------- csgraph : scipy.sparse.csr_matrix @@ -70,7 +74,7 @@ class DiGraph(object): """ - def __init__(self, adj_matrix, weighted=False): + def __init__(self, adj_matrix, weighted=False, node_labels=None): if weighted: dtype = None else: @@ -83,6 +87,10 @@ def __init__(self, adj_matrix, weighted=False): self.n = n # Number of nodes + self._node_labels = None + if node_labels is not None: + self.node_labels = node_labels + self._num_scc = None self._scc_proj = None self._sink_scc_labels = None @@ -95,6 +103,21 @@ def __repr__(self): def __str__(self): return "Directed Graph:\n - n(number of nodes): {n}".format(n=self.n) + @property + def node_labels(self): + return self._node_labels + + @node_labels.setter + def node_labels(self, values): + if len(values) != self.n: + raise ValueError('node_labels must be of length n') + self._node_labels = np.asarray(values) + + def label_nodes(self, list_of_components): + if self.node_labels is not None: + return [self.node_labels[c] for c in list_of_components] + return list_of_components + def _find_scc(self): """ Set ``self._num_scc`` and ``self._scc_proj`` @@ -169,22 +192,42 @@ def sink_scc_labels(self): def num_sink_strongly_connected_components(self): return len(self.sink_scc_labels) - @property - def strongly_connected_components(self): + # strongly_connected_components + def _get_strongly_connected_components(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in range(self.num_strongly_connected_components)] + def get_strongly_connected_components(self, return_labels=True): + if return_labels: + return self.label_nodes(self._get_strongly_connected_components()) + return self._get_strongly_connected_components() + @property - def sink_strongly_connected_components(self): + def strongly_connected_components(self): + return self.get_strongly_connected_components() + + # sink_strongly_connected_components + def _get_sink_strongly_connected_components(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in self.sink_scc_labels.tolist()] + def get_sink_strongly_connected_components(self, return_labels=True): + if return_labels: + return self.label_nodes( + self._get_sink_strongly_connected_components() + ) + return self._get_sink_strongly_connected_components() + + @property + def sink_strongly_connected_components(self): + return self.get_sink_strongly_connected_components() + def _compute_period(self): """ Set ``self._period`` and ``self._cyclic_components_proj``. @@ -255,14 +298,23 @@ def period(self): def is_aperiodic(self): return (self.period == 1) - @property - def cyclic_components(self): + # cyclic_components + def _get_cyclic_components(self): if self.is_aperiodic: return [np.arange(self.n)] else: return [np.where(self._cyclic_components_proj == k)[0] for k in range(self.period)] + def get_cyclic_components(self, return_labels=True): + if return_labels: + return self.label_nodes(self._get_cyclic_components()) + return self._get_cyclic_components() + + @property + def cyclic_components(self): + return self.get_cyclic_components() + def subgraph(self, nodes): """ Return the subgraph consisting of the given nodes and edges @@ -271,7 +323,7 @@ def subgraph(self, nodes): Parameters ---------- nodes : array_like(int, ndim=1) - Array of nodes. + Array of node indices. Returns ------- @@ -282,7 +334,13 @@ def subgraph(self, nodes): adj_matrix = self.csgraph[nodes, :][:, nodes] weighted = True # To copy the dtype - return DiGraph(adj_matrix, weighted=weighted) + + if self.node_labels is not None: + node_labels = self.node_labels[nodes] + else: + node_labels = None + + return DiGraph(adj_matrix, weighted=weighted, node_labels=node_labels) def _csr_matrix_indices(S): From b43039fc48c0bcbe4ed3c4716532b9ef5338d56a Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 15 Mar 2016 21:00:34 -0400 Subject: [PATCH 02/23] Edit and add docstrings --- quantecon/graph_tools.py | 41 +++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index 1fd68a9e0..206939879 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -42,16 +42,18 @@ class DiGraph(object): num_strongly_connected_components : int The number of the strongly connected components. - strongly_connected_components : list(ndarray(int)) + strongly_connected_components : list(ndarray) List of numpy arrays containing the strongly connected - components. + components. Equivalent to calling the method + `get_strongly_connected_components`. num_sink_strongly_connected_components : int The number of the sink strongly connected components. - sink_strongly_connected_components : list(ndarray(int)) + sink_strongly_connected_components : list(ndarray) List of numpy arrays containing the sink strongly connected - components. + components. Equivalent to calling the method + `get_sink_strongly_connected_components`. is_aperiodic : bool Indicate whether the digraph is aperiodic. @@ -60,8 +62,9 @@ class DiGraph(object): The period of the digraph. Defined only for a strongly connected digraph. - cyclic_components : list(ndarray(int)) + cyclic_components : list(ndarray) List of numpy arrays containing the cyclic components. + Equivalent to calling the method `get_cyclic_components`. References ---------- @@ -343,6 +346,34 @@ def subgraph(self, nodes): return DiGraph(adj_matrix, weighted=weighted, node_labels=node_labels) +_get_method_docstr = \ +""" +Return a list of numpy arrays containing the {components}. + +Parameters +---------- +return_labels : bool(optional, default=True) + Whether to annotate the returned nodes with `node_labels`. + +Returns +------- +list(ndarray) + If `return_labels=True`, and if `node_labels` is not None, + each ndarray contains the node labels, and the node indices + (integers) otherwise. + +""" + +DiGraph.get_strongly_connected_components.__doc__ = \ + _get_method_docstr.format(components='strongly connected components') + +DiGraph.get_sink_strongly_connected_components.__doc__ = \ + _get_method_docstr.format(components='sink strongly connected components') + +DiGraph.get_cyclic_components.__doc__ = \ + _get_method_docstr.format(components='cyclic components') + + def _csr_matrix_indices(S): """ Generate the indices of nonzero entries of a csr_matrix S From 0176d3e4360b820141299ef14a2339264afe82e4 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 15 Mar 2016 22:09:00 -0400 Subject: [PATCH 03/23] Add tests for node labels --- quantecon/tests/test_graph_tools.py | 54 +++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/quantecon/tests/test_graph_tools.py b/quantecon/tests/test_graph_tools.py index dc83afb0c..b7d766ae7 100644 --- a/quantecon/tests/test_graph_tools.py +++ b/quantecon/tests/test_graph_tools.py @@ -228,6 +228,60 @@ def test_subgraph_weighted(): ) +def test_node_labels_connected_components(): + adj_matrix = [[1, 0, 0], [1, 0, 0], [0, 0, 1]] + node_labels = np.array(['a', 'b', 'c']) + g = DiGraph(adj_matrix, node_labels=node_labels) + + sccs = [[0], [1], [2]] + sink_sccs = [[0], [2]] + + methods = ['get_strongly_connected_components', + 'get_sink_strongly_connected_components'] + for method, components_ind in zip(methods, [sccs, sink_sccs]): + for return_labels in [True, False]: + if return_labels: + components = [node_labels[i] for i in components_ind] + else: + components = components_ind + list_of_array_equal( + sorted(getattr(g, method)(return_labels), key=lambda x: x[0]), + sorted(components, key=lambda x: x[0]) + ) + + +def test_node_labels_cyclic_components(): + adj_matrix = [[0, 1], [1, 0]] + node_labels = np.array(['a', 'b']) + g = DiGraph(adj_matrix, node_labels=node_labels) + + cyclic_components = [[0], [1]] + + methods = ['get_cyclic_components'] + for method, components_ind in zip(methods, [cyclic_components]): + for return_labels in [True, False]: + if return_labels: + components = [node_labels[i] for i in components_ind] + else: + components = components_ind + list_of_array_equal( + sorted(getattr(g, method)(return_labels), key=lambda x: x[0]), + sorted(components, key=lambda x: x[0]) + ) + + +def test_node_labels_subgraph(): + adj_matrix = [[0, 1, 0], [0, 0, 1], [1, 0, 0]] + node_labels = np.array(['a', 'b', 'c']) + g = DiGraph(adj_matrix, node_labels=node_labels) + nodes = [1, 2] + + assert_array_equal( + g.subgraph(nodes).node_labels, + node_labels[nodes] + ) + + @raises(ValueError) def test_raises_value_error_non_sym(): """Test with non symmetric input""" From 3c103d68f20f79b2ac244cc8d202fa5faeba1b03 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 12:01:39 -0400 Subject: [PATCH 04/23] Add decorator `annotate_nodes` to get_*_components methods Properties return indices as before --- quantecon/graph_tools.py | 75 ++++++++++++++++++---------------------- 1 file changed, 33 insertions(+), 42 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index 206939879..d71887e83 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -12,6 +12,16 @@ from fractions import gcd +# Decorator for get_*_components methods +def annotate_nodes(func): + def new_func(self, return_labels=True): + list_of_components, return_labels = func(self, return_labels) + if return_labels and (self.node_labels is not None): + return [self.node_labels[c] for c in list_of_components] + return list_of_components + return new_func + + class DiGraph(object): r""" Class for a directed graph. It stores useful information about the @@ -42,18 +52,16 @@ class DiGraph(object): num_strongly_connected_components : int The number of the strongly connected components. - strongly_connected_components : list(ndarray) - List of numpy arrays containing the strongly connected - components. Equivalent to calling the method - `get_strongly_connected_components`. + strongly_connected_components : list(ndarray(int)) + List of numpy arrays containing the indices of the strongly + connected components. num_sink_strongly_connected_components : int The number of the sink strongly connected components. - sink_strongly_connected_components : list(ndarray) - List of numpy arrays containing the sink strongly connected - components. Equivalent to calling the method - `get_sink_strongly_connected_components`. + sink_strongly_connected_components : list(ndarray(int)) + List of numpy arrays containing the indices of the sink strongly + connected components. is_aperiodic : bool Indicate whether the digraph is aperiodic. @@ -62,9 +70,9 @@ class DiGraph(object): The period of the digraph. Defined only for a strongly connected digraph. - cyclic_components : list(ndarray) - List of numpy arrays containing the cyclic components. - Equivalent to calling the method `get_cyclic_components`. + cyclic_components : list(ndarray(int)) + List of numpy arrays containing the indices of the cyclic + components. References ---------- @@ -195,41 +203,29 @@ def sink_scc_labels(self): def num_sink_strongly_connected_components(self): return len(self.sink_scc_labels) - # strongly_connected_components - def _get_strongly_connected_components(self): + @property + def strongly_connected_components(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in range(self.num_strongly_connected_components)] - def get_strongly_connected_components(self, return_labels=True): - if return_labels: - return self.label_nodes(self._get_strongly_connected_components()) - return self._get_strongly_connected_components() + @annotate_nodes + def get_strongly_connected_components(self, return_labels): + return self.strongly_connected_components, return_labels @property - def strongly_connected_components(self): - return self.get_strongly_connected_components() - - # sink_strongly_connected_components - def _get_sink_strongly_connected_components(self): + def sink_strongly_connected_components(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in self.sink_scc_labels.tolist()] - def get_sink_strongly_connected_components(self, return_labels=True): - if return_labels: - return self.label_nodes( - self._get_sink_strongly_connected_components() - ) - return self._get_sink_strongly_connected_components() - - @property - def sink_strongly_connected_components(self): - return self.get_sink_strongly_connected_components() + @annotate_nodes + def get_sink_strongly_connected_components(self, return_labels): + return self.sink_strongly_connected_components, return_labels def _compute_period(self): """ @@ -301,22 +297,17 @@ def period(self): def is_aperiodic(self): return (self.period == 1) - # cyclic_components - def _get_cyclic_components(self): + @property + def cyclic_components(self): if self.is_aperiodic: return [np.arange(self.n)] else: return [np.where(self._cyclic_components_proj == k)[0] for k in range(self.period)] - def get_cyclic_components(self, return_labels=True): - if return_labels: - return self.label_nodes(self._get_cyclic_components()) - return self._get_cyclic_components() - - @property - def cyclic_components(self): - return self.get_cyclic_components() + @annotate_nodes + def get_cyclic_components(self, return_labels): + return self.cyclic_components, return_labels def subgraph(self, nodes): """ From 01fb48bced9f1dc1dc2c8db9badfd5086d152b43 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 12:11:26 -0400 Subject: [PATCH 05/23] Delete unused method --- quantecon/graph_tools.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index d71887e83..b6d7b484f 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -124,11 +124,6 @@ def node_labels(self, values): raise ValueError('node_labels must be of length n') self._node_labels = np.asarray(values) - def label_nodes(self, list_of_components): - if self.node_labels is not None: - return [self.node_labels[c] for c in list_of_components] - return list_of_components - def _find_scc(self): """ Set ``self._num_scc`` and ``self._scc_proj`` From 3d243f19bad1e4bfa2e88da7b0803b9187b211b3 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 13:53:52 -0400 Subject: [PATCH 06/23] Allow setting node_labels to None --- quantecon/graph_tools.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index b6d7b484f..be348bd73 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -98,9 +98,8 @@ def __init__(self, adj_matrix, weighted=False, node_labels=None): self.n = n # Number of nodes - self._node_labels = None - if node_labels is not None: - self.node_labels = node_labels + # Call the setter method + self.node_labels = node_labels self._num_scc = None self._scc_proj = None @@ -120,9 +119,12 @@ def node_labels(self): @node_labels.setter def node_labels(self, values): - if len(values) != self.n: - raise ValueError('node_labels must be of length n') - self._node_labels = np.asarray(values) + if values is None: + self._node_labels = None + else: + if len(values) != self.n: + raise ValueError('node_labels must be of length n') + self._node_labels = np.asarray(values) def _find_scc(self): """ From 43632eb537dd5c2770cb5b1cb4074b80fc9ae52a Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 13:54:51 -0400 Subject: [PATCH 07/23] Add state values to MarkovChain --- quantecon/markov/core.py | 87 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 82 insertions(+), 5 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index 4f887c7f0..2df3be1aa 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -108,6 +108,11 @@ class MarkovChain(object): P : array_like or scipy sparse matrix (float, ndim=2) The transition matrix. Must be of shape n x n. + state_values : array_like(ndim=1, default=None) + Array_like of length n containing the values associated with + each states. If None, the values default to integers 0 through + n-1. + Attributes ---------- P : ndarray or scipy.sparse.csr_matrix (float, ndim=2) @@ -149,7 +154,7 @@ class MarkovChain(object): """ - def __init__(self, P): + def __init__(self, P, state_values=None): if sparse.issparse(P): # Sparse matrix self.P = sparse.csr_matrix(P) self.is_sparse = True @@ -180,6 +185,9 @@ def __init__(self, P): if not np.allclose(row_sums, np.ones(self.n)): raise ValueError('The rows of P must sum to 1') + # Call the setter method + self._state_values = state_values + # To analyze the structure of P as a directed graph self._digraph = None @@ -199,10 +207,23 @@ def __repr__(self): def __str__(self): return str(self.__repr__) + @property + def state_values(self): + return self._state_values + + @state_values.setter + def state_values(self, values): + if values is None: + self._state_values = None + else: + if len(values) != self.n: + raise ValueError('state_values must be of length n') + self._state_values = np.asarray(values) + @property def digraph(self): if self._digraph is None: - self._digraph = DiGraph(self.P) + self._digraph = DiGraph(self.P, node_labels=self.state_values) return self._digraph @property @@ -217,6 +238,11 @@ def num_communication_classes(self): def communication_classes(self): return self.digraph.strongly_connected_components + def get_communication_classes(self, return_values): + return self.digraph.get_strongly_connected_components( + return_labels=return_values + ) + @property def num_recurrent_classes(self): return self.digraph.num_sink_strongly_connected_components @@ -225,6 +251,11 @@ def num_recurrent_classes(self): def recurrent_classes(self): return self.digraph.sink_strongly_connected_components + def get_recurrent_classes(self, return_values): + return self.digraph.sink_strongly_connected_components( + return_labels=return_values + ) + @property def is_aperiodic(self): if self.is_irreducible: @@ -256,6 +287,14 @@ def cyclic_classes(self): else: return self.digraph.cyclic_components + def get_cyclic_classes(self, return_values): + if not self.is_irreducible: + raise NotImplementedError( + 'Not defined for a reducible Markov chain' + ) + else: + return self.digraph.cyclic_components(return_labels=return_values) + def _compute_stationary(self): """ Store the stationary distributions in self._stationary_distributions. @@ -311,7 +350,8 @@ def cdfs1d(self): self._cdfs1d = cdfs1d return self._cdfs1d - def simulate(self, ts_length, init=None, num_reps=None, random_state=None): + def simulate(self, ts_length, init=None, num_reps=None, return_values=True, + random_state=None): """ Simulate time series of state transitions. @@ -328,6 +368,9 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None): num_reps : scalar(int), optional(default=None) Number of repetitions of simulation. + return_values : bool(optional, default=True) + Whether to annotate the returned states with state_values. + random_state : scalar(int) or np.random.RandomState, optional(default=None) Random seed (integer) or np.random.RandomState instance to @@ -337,13 +380,15 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None): Returns ------- - X : ndarray(int, ndim=1 or 2) + X : ndarray(ndim=1 or 2) Array containing the sample path(s), of shape (ts_length,) if init is a scalar (integer) or None and num_reps is None; of shape (k, ts_length) otherwise, where k = len(init) if (init, num_reps) = (array, None), k = num_reps if (init, num_reps) = (int or None, int), and k = len(init)*num_reps - if (init, num_reps) = (array, int). + if (init, num_reps) = (array, int). If return_values=True, + and if state_values is not None, the array elements are the + state values, and the state indices (integers) otherwise. """ random_state = check_random_state(random_state) @@ -398,6 +443,10 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None): random_values, out=X ) + # Annotate states + if self.state_values is not None: + X = self.state_values[X] + if dim == 1: return X[0] else: @@ -490,6 +539,34 @@ def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states, jit(nopython=True)(_generate_sample_paths_sparse) +_get_method_docstr = \ +""" +Return a list of numpy arrays containing the {classes}. + +Parameters +---------- +return_values : bool(optional, default=True) + Whether to annotate the returned states with `state_values`. + +Returns +------- +list(ndarray) + If `return_values=True`, and if `state_values` is not None, + each ndarray contains the state values, and the state indices + (integers) otherwise. + +""" + +MarkovChain.get_communication_classes.__doc__ = \ + _get_method_docstr.format(classes='communication classes') + +MarkovChain.get_recurrent_classes.__doc__ = \ + _get_method_docstr.format(classes='recurrent classes') + +MarkovChain.get_cyclic_classes.__doc__ = \ + _get_method_docstr.format(classes='cyclic classes') + + def mc_compute_stationary(P): """ Computes stationary distributions of P, one for each recurrent From b0ad892f863d2b122b4ad5fc487da1307fa2c837 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 16:38:45 -0400 Subject: [PATCH 08/23] Fix bugs --- quantecon/markov/core.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index 2df3be1aa..0c8b54a7e 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -252,7 +252,7 @@ def recurrent_classes(self): return self.digraph.sink_strongly_connected_components def get_recurrent_classes(self, return_values): - return self.digraph.sink_strongly_connected_components( + return self.digraph.get_sink_strongly_connected_components( return_labels=return_values ) @@ -293,7 +293,9 @@ def get_cyclic_classes(self, return_values): 'Not defined for a reducible Markov chain' ) else: - return self.digraph.cyclic_components(return_labels=return_values) + return self.digraph.get_cyclic_components( + return_labels=return_values + ) def _compute_stationary(self): """ @@ -444,7 +446,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, ) # Annotate states - if self.state_values is not None: + if return_values and (self.state_values is not None): X = self.state_values[X] if dim == 1: From b5cacb9d4487ec50953806508ae580346c887249 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Thu, 17 Mar 2016 17:00:48 -0400 Subject: [PATCH 09/23] Add tests for state values --- quantecon/markov/tests/test_core.py | 79 ++++++++++++++++++++++++++++- 1 file changed, 78 insertions(+), 1 deletion(-) diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 39ba0e653..6b31dbfcb 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -1,5 +1,5 @@ """ -Tests for mc_tools.py +Tests for markov/core.py Functions --------- @@ -381,6 +381,83 @@ def test_mc_sample_path_lln(): ok_(np.abs(frequency_1 - stationary_dist[1]) < tol) +class TestMCStateValues: + def setUp(self): + self.state_values = np.array([[0, 1], [2, 3], [4, 5]]) + + self.mc_reducible_dict = { + 'mc': MarkovChain([[1, 0, 0], [1, 0, 0], [0, 0, 1]], + state_values=self.state_values), + 'coms': [[0], [1], [2]], + 'recs': [[0], [2]] + } + + self.mc_periodic_dict = { + 'mc': MarkovChain([[0, 1, 0], [0, 0, 1], [1, 0, 0]], + state_values=self.state_values), + 'cycs': [[0], [1], [2]] + } + + def test_com_rec_classes(self): + mc = self.mc_reducible_dict['mc'] + coms = self.mc_reducible_dict['coms'] + recs = self.mc_reducible_dict['recs'] + methods = ['get_communication_classes', + 'get_recurrent_classes'] + for method, classes_ind in zip(methods, [coms, recs]): + for return_values in [True, False]: + if return_values: + classes = [self.state_values[i] for i in classes_ind] + key = lambda x: x[0, 0] + else: + classes = classes_ind + key = lambda x: x[0] + list_of_array_equal( + sorted(getattr(mc, method)(return_values), key=key), + sorted(classes, key=key) + ) + + def test_cyc_classes(self): + mc = self.mc_periodic_dict['mc'] + cycs = self.mc_periodic_dict['cycs'] + methods = ['get_cyclic_classes'] + for method, classes_ind in zip(methods, [cycs]): + for return_values in [True, False]: + if return_values: + classes = [self.state_values[i] for i in classes_ind] + key = lambda x: x[0, 0] + else: + classes = classes_ind + key = lambda x: x[0] + list_of_array_equal( + sorted(getattr(mc, method)(return_values), key=key), + sorted(classes, key=key) + ) + + def test_simulate(self): + # Deterministic mc + mc = self.mc_periodic_dict['mc'] + ts_length = 6 + + init = 0 + for return_values in [True, False]: + X = mc.simulate(ts_length, init=init, return_values=return_values) + X_expected = np.arange(init, init+ts_length)%mc.n + if return_values: + X_expected = self.state_values[X_expected] + assert_array_equal(X, X_expected) + + inits = [1, 2] + for return_values in [True, False]: + X = mc.simulate(ts_length, init=inits, return_values=return_values) + X_expected = np.array( + [np.arange(init, init+ts_length)%mc.n for init in inits] + ) + if return_values: + X_expected = self.state_values[X_expected] + assert_array_equal(X, X_expected) + + @raises(ValueError) def test_raises_value_error_non_2dim(): """Test with non 2dim input""" From 84fdbaafedfe89061fcde8fe60be206716435357 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Fri, 18 Mar 2016 16:56:50 +0900 Subject: [PATCH 10/23] Add a testcase --- quantecon/markov/tests/test_core.py | 37 ++++++++++++++++------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 6b31dbfcb..3c10216f8 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -395,27 +395,30 @@ def setUp(self): self.mc_periodic_dict = { 'mc': MarkovChain([[0, 1, 0], [0, 0, 1], [1, 0, 0]], state_values=self.state_values), + 'coms': [[0, 1, 2]], + 'recs': [[0, 1, 2]], 'cycs': [[0], [1], [2]] } def test_com_rec_classes(self): - mc = self.mc_reducible_dict['mc'] - coms = self.mc_reducible_dict['coms'] - recs = self.mc_reducible_dict['recs'] - methods = ['get_communication_classes', - 'get_recurrent_classes'] - for method, classes_ind in zip(methods, [coms, recs]): - for return_values in [True, False]: - if return_values: - classes = [self.state_values[i] for i in classes_ind] - key = lambda x: x[0, 0] - else: - classes = classes_ind - key = lambda x: x[0] - list_of_array_equal( - sorted(getattr(mc, method)(return_values), key=key), - sorted(classes, key=key) - ) + for mc_dict in [self.mc_reducible_dict, self.mc_periodic_dict]: + mc = mc_dict['mc'] + coms = mc_dict['coms'] + recs = mc_dict['recs'] + methods = ['get_communication_classes', + 'get_recurrent_classes'] + for method, classes_ind in zip(methods, [coms, recs]): + for return_values in [True, False]: + if return_values: + classes = [self.state_values[i] for i in classes_ind] + key = lambda x: x[0, 0] + else: + classes = classes_ind + key = lambda x: x[0] + list_of_array_equal( + sorted(getattr(mc, method)(return_values), key=key), + sorted(classes, key=key) + ) def test_cyc_classes(self): mc = self.mc_periodic_dict['mc'] From 457578b1013417a2839d61d8e9e9bfa2a9e9c122 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sat, 19 Mar 2016 20:53:32 +0900 Subject: [PATCH 11/23] MarkovChain: Remove util.numba_installed --- quantecon/markov/core.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index 0c8b54a7e..b8854ac43 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -86,14 +86,14 @@ """ from __future__ import division import numbers +from fractions import gcd import numpy as np from scipy import sparse -from fractions import gcd +from numba import jit + from .gth_solve import gth_solve from ..graph_tools import DiGraph - -# -Check if Numba is Available- # -from ..util import searchsorted, check_random_state, numba_installed, jit +from ..util import searchsorted, check_random_state class MarkovChain(object): @@ -455,6 +455,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, return X +@jit(nopython=True) def _generate_sample_paths(P_cdfs, init_states, random_values, out): """ Generate num_reps sample paths of length ts_length, where num_reps = @@ -478,7 +479,7 @@ def _generate_sample_paths(P_cdfs, init_states, random_values, out): Notes ----- - This routine is jit-complied if the module Numba is vailable. + This routine is jit-complied by Numba. """ num_reps, ts_length = out.shape @@ -488,10 +489,8 @@ def _generate_sample_paths(P_cdfs, init_states, random_values, out): for t in range(ts_length-1): out[i, t+1] = searchsorted(P_cdfs[out[i, t]], random_values[i, t]) -if numba_installed: - _generate_sample_paths = jit(nopython=True)(_generate_sample_paths) - +@jit(nopython=True) def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states, random_values, out): """ @@ -524,7 +523,7 @@ def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states, Notes ----- - This routine is jit-complied if the module Numba is vailable. + This routine is jit-complied by Numba. """ num_reps, ts_length = out.shape @@ -536,10 +535,6 @@ def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states, random_values[i, t]) out[i, t+1] = indices[indptr[out[i, t]]+k] -if numba_installed: - _generate_sample_paths_sparse = \ - jit(nopython=True)(_generate_sample_paths_sparse) - _get_method_docstr = \ """ From 8850b31a72f7fa478425be3b4bb0762ba9090cdf Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sat, 19 Mar 2016 20:54:18 +0900 Subject: [PATCH 12/23] MarkovChain: Fix bug --- quantecon/markov/core.py | 2 +- quantecon/markov/tests/test_core.py | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index b8854ac43..e9034ecbc 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -186,7 +186,7 @@ def __init__(self, P, state_values=None): raise ValueError('The rows of P must sum to 1') # Call the setter method - self._state_values = state_values + self.state_values = state_values # To analyze the structure of P as a directed graph self._digraph = None diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 3c10216f8..0f5679241 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -383,18 +383,19 @@ def test_mc_sample_path_lln(): class TestMCStateValues: def setUp(self): - self.state_values = np.array([[0, 1], [2, 3], [4, 5]]) + state_values = [[0, 1], [2, 3], [4, 5]] # Pass python list + self.state_values = np.array(state_values) self.mc_reducible_dict = { 'mc': MarkovChain([[1, 0, 0], [1, 0, 0], [0, 0, 1]], - state_values=self.state_values), + state_values=state_values), 'coms': [[0], [1], [2]], 'recs': [[0], [2]] } self.mc_periodic_dict = { 'mc': MarkovChain([[0, 1, 0], [0, 0, 1], [1, 0, 0]], - state_values=self.state_values), + state_values=state_values), 'coms': [[0, 1, 2]], 'recs': [[0, 1, 2]], 'cycs': [[0], [1], [2]] From 8d3176ccebb57b95407fcfd8d53939941919d951 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 21 Mar 2016 00:00:39 +0900 Subject: [PATCH 13/23] MarkovChain: Set default `return_values=True` --- quantecon/markov/core.py | 6 +++--- quantecon/markov/tests/test_core.py | 16 ++++++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index e9034ecbc..f6812449a 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -238,7 +238,7 @@ def num_communication_classes(self): def communication_classes(self): return self.digraph.strongly_connected_components - def get_communication_classes(self, return_values): + def get_communication_classes(self, return_values=True): return self.digraph.get_strongly_connected_components( return_labels=return_values ) @@ -251,7 +251,7 @@ def num_recurrent_classes(self): def recurrent_classes(self): return self.digraph.sink_strongly_connected_components - def get_recurrent_classes(self, return_values): + def get_recurrent_classes(self, return_values=True): return self.digraph.get_sink_strongly_connected_components( return_labels=return_values ) @@ -287,7 +287,7 @@ def cyclic_classes(self): else: return self.digraph.cyclic_components - def get_cyclic_classes(self, return_values): + def get_cyclic_classes(self, return_values=True): if not self.is_irreducible: raise NotImplementedError( 'Not defined for a reducible Markov chain' diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 0f5679241..24507a291 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -420,6 +420,14 @@ def test_com_rec_classes(self): sorted(getattr(mc, method)(return_values), key=key), sorted(classes, key=key) ) + # Check that the default of return_values is True + classes = [self.state_values[i] for i in classes_ind] + key = lambda x: x[0, 0] + list_of_array_equal( + sorted(getattr(mc, method)(), key=key), + sorted(classes, key=key) + ) + def test_cyc_classes(self): mc = self.mc_periodic_dict['mc'] @@ -437,6 +445,14 @@ def test_cyc_classes(self): sorted(getattr(mc, method)(return_values), key=key), sorted(classes, key=key) ) + # Check that the default of return_values is True + classes = [self.state_values[i] for i in classes_ind] + key = lambda x: x[0, 0] + list_of_array_equal( + sorted(getattr(mc, method)(), key=key), + sorted(classes, key=key) + ) + def test_simulate(self): # Deterministic mc From ecfbc8f65f25468e0e61f0fafdbf553f8a9573ca Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 28 Mar 2016 22:14:40 +0900 Subject: [PATCH 14/23] Disallow non-homogeneous node_labels/state_values --- quantecon/graph_tools.py | 20 ++++++++++++++------ quantecon/markov/core.py | 21 ++++++++++++++------- quantecon/markov/tests/test_core.py | 6 ++++++ quantecon/tests/test_graph_tools.py | 8 +++++++- 4 files changed, 41 insertions(+), 14 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index be348bd73..0e966729c 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -37,9 +37,10 @@ class DiGraph(object): weighted : bool, optional(default=False) Whether to treat `adj_matrix` as a weighted adjacency matrix. - node_labels : array_like(ndim=1, default=None) - Array_like of length n containing the label associated with each - node. If None, the labels default to integers 0 through n-1. + node_labels : array_like(default=None) + Array_like of length n containing the labels associated with the + states, which must be homogeneous in type. If None, the labels + default to integers 0 through n-1. Attributes ---------- @@ -122,9 +123,16 @@ def node_labels(self, values): if values is None: self._node_labels = None else: - if len(values) != self.n: - raise ValueError('node_labels must be of length n') - self._node_labels = np.asarray(values) + values = np.asarray(values) + if (values.ndim < 1) or (values.shape[0] != self.n): + raise ValueError( + 'node_labels must be an array_like of length n' + ) + if np.issubdtype(values.dtype, np.object_): + raise ValueError( + 'data in node_labels must be homogeneous in type' + ) + self._node_labels = values def _find_scc(self): """ diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index f6812449a..b4c513ddc 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -108,10 +108,10 @@ class MarkovChain(object): P : array_like or scipy sparse matrix (float, ndim=2) The transition matrix. Must be of shape n x n. - state_values : array_like(ndim=1, default=None) - Array_like of length n containing the values associated with - each states. If None, the values default to integers 0 through - n-1. + state_values : array_like(default=None) + Array_like of length n containing the values associated with the + states, which must be homogeneous in type. If None, the values + default to integers 0 through n-1. Attributes ---------- @@ -216,9 +216,16 @@ def state_values(self, values): if values is None: self._state_values = None else: - if len(values) != self.n: - raise ValueError('state_values must be of length n') - self._state_values = np.asarray(values) + values = np.asarray(values) + if (values.ndim < 1) or (values.shape[0] != self.n): + raise ValueError( + 'state_values must be an array_like of length n' + ) + if np.issubdtype(values.dtype, np.object_): + raise ValueError( + 'data in state_values must be homogeneous in type' + ) + self._state_values = values @property def digraph(self): diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 24507a291..c5b31fb9b 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -517,6 +517,12 @@ def test_raises_value_error_simulate_init_out_of_range(): assert_raises(ValueError, mc.simulate, ts_length, init=[0, -(n+1)]) +def test_raises_non_homogeneous_state_values(): + P = [[0.4, 0.6], [0.2, 0.8]] + state_values = [(0, 1), 2] + assert_raises(ValueError, MarkovChain, P, state_values=state_values) + + if __name__ == '__main__': import sys import nose diff --git a/quantecon/tests/test_graph_tools.py b/quantecon/tests/test_graph_tools.py index b7d766ae7..f646e5c08 100644 --- a/quantecon/tests/test_graph_tools.py +++ b/quantecon/tests/test_graph_tools.py @@ -7,7 +7,7 @@ """ import sys import numpy as np -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_raises import nose from nose.tools import eq_, ok_, raises @@ -288,6 +288,12 @@ def test_raises_value_error_non_sym(): g = DiGraph(np.array([[0.4, 0.6]])) +def test_raises_non_homogeneous_node_labels(): + adj_matrix = [[1, 0], [0, 1]] + node_labels = [(0, 1), 2] + assert_raises(ValueError, DiGraph, adj_matrix, node_labels=node_labels) + + if __name__ == '__main__': argv = sys.argv[:] argv.append('--verbose') From 5bc78d5706462cde82dcba09dfec114c2a3e4680 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 28 Mar 2016 22:38:45 +0900 Subject: [PATCH 15/23] MarkovChain: Add get_index --- quantecon/markov/core.py | 38 +++++++++++++++++++++++++++++ quantecon/markov/tests/test_core.py | 20 +++++++++++++-- 2 files changed, 56 insertions(+), 2 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index b4c513ddc..c25cc3cac 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -227,6 +227,44 @@ def state_values(self, values): ) self._state_values = values + def get_index(self, value): + """ + Return the index of the given value in state_values. + + Parameters + ---------- + value + Value to get the index for. + + Returns + ------- + idx : int + Index of the value. + + """ + error_msg = 'value {0} not found'.format(repr(value)) + + if self.state_values is None: + if isinstance(value, numbers.Integral) and (0 <= value < self.n): + return value + else: + raise ValueError(error_msg) + + # if self.state_values is not None: + if self.state_values.ndim == 1: + try: + idx = np.where(self.state_values == value)[0][0] + return idx + except IndexError: + raise ValueError(error_msg) + else: + idx = 0 + while idx < self.n: + if np.array_equal(self.state_values[idx], value): + return idx + idx += 1 + raise ValueError(error_msg) + @property def digraph(self): if self._digraph is None: diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index c5b31fb9b..512e3c4b8 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -428,7 +428,6 @@ def test_com_rec_classes(self): sorted(classes, key=key) ) - def test_cyc_classes(self): mc = self.mc_periodic_dict['mc'] cycs = self.mc_periodic_dict['cycs'] @@ -453,7 +452,6 @@ def test_cyc_classes(self): sorted(classes, key=key) ) - def test_simulate(self): # Deterministic mc mc = self.mc_periodic_dict['mc'] @@ -478,6 +476,24 @@ def test_simulate(self): assert_array_equal(X, X_expected) +def test_get_index(): + P = [[0.4, 0.6], [0.2, 0.8]] + mc = MarkovChain(P) + + eq_(mc.get_index(0), 0) + eq_(mc.get_index(1), 1) + assert_raises(ValueError, mc.get_index, 2) + + mc.state_values = [1, 2] + eq_(mc.get_index(1), 0) + eq_(mc.get_index(2), 1) + assert_raises(ValueError, mc.get_index, 0) + + mc.state_values = [[1, 2], [3, 4]] + eq_(mc.get_index([1, 2]), 0) + assert_raises(ValueError, mc.get_index, 1) + + @raises(ValueError) def test_raises_value_error_non_2dim(): """Test with non 2dim input""" From bfc91bb69443e2e234643ecc26fb2c05fe24f320 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 29 Mar 2016 00:20:05 +0900 Subject: [PATCH 16/23] MarkovChain.get_index: Accept array_like of values --- quantecon/markov/core.py | 42 ++++++++++++++++++++++++++--- quantecon/markov/tests/test_core.py | 5 ++++ 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index c25cc3cac..2a0d49d37 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -229,7 +229,43 @@ def state_values(self, values): def get_index(self, value): """ - Return the index of the given value in state_values. + Return the index (or indices) of the given value (or values) in + `state_values`. + + Parameters + ---------- + value + Value(s) to get the index (indices) for. + + Returns + ------- + idx : int or ndarray(int) + Index of `value` if `value` is a single state value; array + of indices if `value` is an array_like of state values. + + """ + if self.state_values is None: + state_values_ndim = 1 + else: + state_values_ndim = self.state_values.ndim + + values = np.asarray(value) + + if values.ndim <= state_values_ndim - 1: + return self._get_index(value) + elif values.ndim == state_values_ndim: # array of values + k = values.shape[0] + idx = np.empty(k, dtype=int) + for i in range(k): + idx[i] = self._get_index(values[i]) + return idx + else: + raise ValueError('invalid value') + + + def _get_index(self, value): + """ + Return the index of the given value in `state_values`. Parameters ---------- @@ -239,10 +275,10 @@ def get_index(self, value): Returns ------- idx : int - Index of the value. + Index of `value`. """ - error_msg = 'value {0} not found'.format(repr(value)) + error_msg = 'value {0} not found'.format(value) if self.state_values is None: if isinstance(value, numbers.Integral) and (0 <= value < self.n): diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 512e3c4b8..44410b593 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -483,15 +483,20 @@ def test_get_index(): eq_(mc.get_index(0), 0) eq_(mc.get_index(1), 1) assert_raises(ValueError, mc.get_index, 2) + assert_array_equal(mc.get_index([1, 0]), [1, 0]) + assert_raises(ValueError, mc.get_index, [[1]]) mc.state_values = [1, 2] eq_(mc.get_index(1), 0) eq_(mc.get_index(2), 1) assert_raises(ValueError, mc.get_index, 0) + assert_array_equal(mc.get_index([2, 1]), [1, 0]) + assert_raises(ValueError, mc.get_index, [[1]]) mc.state_values = [[1, 2], [3, 4]] eq_(mc.get_index([1, 2]), 0) assert_raises(ValueError, mc.get_index, 1) + assert_array_equal(mc.get_index([[3, 4], [1, 2]]), [1, 0]) @raises(ValueError) From 193f625911010eb71d93d7749ffa8dc67cd0326b Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 29 Mar 2016 11:05:07 +0900 Subject: [PATCH 17/23] Fix typo --- quantecon/markov/core.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index 2a0d49d37..ca8f5a761 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -477,7 +477,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, random_state = check_random_state(random_state) dim = 1 # Dimension of the returned array: 1 or 2 - msg_out_or_range = 'index {init} is out of the state space' + msg_out_of_range = 'index {init} is out of the state space' try: k = len(init) # init is an array @@ -488,7 +488,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, idx = np.where( (init_states >= self.n) + (init_states < -self.n) )[0][0] - raise ValueError(msg_out_or_range.format(init=idx)) + raise ValueError(msg_out_of_range.format(init=idx)) if num_reps is not None: k *= num_reps init_states = np.tile(init_states, num_reps) @@ -502,7 +502,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, elif isinstance(init, numbers.Integral): # Check init is in the state space if init >= self.n or init < -self.n: - raise ValueError(msg_out_or_range.format(init=init)) + raise ValueError(msg_out_of_range.format(init=init)) init_states = np.ones(k, dtype=int) * init else: raise ValueError( From d825975486cabe788730e86b52a3740533f136c6 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 29 Mar 2016 11:53:55 +0900 Subject: [PATCH 18/23] MarkovChain: Add `simulate_values` and remove `return_values` from `simulate` --- quantecon/markov/core.py | 56 ++++++++++++++++++++++------- quantecon/markov/tests/test_core.py | 30 ++++++++-------- 2 files changed, 60 insertions(+), 26 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index ca8f5a761..cd035e2c5 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -433,8 +433,7 @@ def cdfs1d(self): self._cdfs1d = cdfs1d return self._cdfs1d - def simulate(self, ts_length, init=None, num_reps=None, return_values=True, - random_state=None): + def simulate(self, ts_length, init=None, num_reps=None, random_state=None): """ Simulate time series of state transitions. @@ -451,9 +450,6 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, num_reps : scalar(int), optional(default=None) Number of repetitions of simulation. - return_values : bool(optional, default=True) - Whether to annotate the returned states with state_values. - random_state : scalar(int) or np.random.RandomState, optional(default=None) Random seed (integer) or np.random.RandomState instance to @@ -469,9 +465,7 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, of shape (k, ts_length) otherwise, where k = len(init) if (init, num_reps) = (array, None), k = num_reps if (init, num_reps) = (int or None, int), and k = len(init)*num_reps - if (init, num_reps) = (array, int). If return_values=True, - and if state_values is not None, the array elements are the - state values, and the state indices (integers) otherwise. + if (init, num_reps) = (array, int). """ random_state = check_random_state(random_state) @@ -526,15 +520,53 @@ def simulate(self, ts_length, init=None, num_reps=None, return_values=True, random_values, out=X ) - # Annotate states - if return_values and (self.state_values is not None): - X = self.state_values[X] - if dim == 1: return X[0] else: return X + def simulate_values(self, ts_length, init_value=None, num_reps=None, + random_state=None): + """ + Simulate time series of state transitions, but return the values + in `state_values` instead of the state indices. + + Parameters + ---------- + ts_length : scalar(int) + Length of each simulation. + + init_value : scalar or array_like, optional(default=None) + Initial state values(s). If None, the initial state is + randomly drawn. + + num_reps : scalar(int), optional(default=None) + Number of repetitions of simulation. + + random_state : scalar(int) or np.random.RandomState, + optional(default=None) + Random seed (integer) or np.random.RandomState instance to + set the initial state of the random number generator for + reproducibility. If None, a randomly initialized RandomState + is used. + + Returns + ------- + X : ndarray(ndim=1 or 2) + Array containing the state values of the sample path(s). See + the `simulate` method for more information. + + """ + init_idx = self.get_index(init_value) + X = self.simulate(ts_length, init=init_idx, num_reps=num_reps, + random_state=random_state) + + # Annotate states + if self.state_values is not None: + X = self.state_values[X] + + return X + @jit(nopython=True) def _generate_sample_paths(P_cdfs, init_states, random_values, out): diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 44410b593..39e1c49d2 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -457,22 +457,24 @@ def test_simulate(self): mc = self.mc_periodic_dict['mc'] ts_length = 6 - init = 0 - for return_values in [True, False]: - X = mc.simulate(ts_length, init=init, return_values=return_values) - X_expected = np.arange(init, init+ts_length)%mc.n - if return_values: - X_expected = self.state_values[X_expected] + methods = ['simulate', 'simulate_values'] + + init_idx = 0 + inits = [init_idx, self.state_values[init_idx]] + path = np.arange(init_idx, init_idx+ts_length)%mc.n + paths = [path, self.state_values[path]] + for method, init, X_expected in zip(methods, inits, paths): + X = getattr(mc, method)(ts_length, init) assert_array_equal(X, X_expected) - inits = [1, 2] - for return_values in [True, False]: - X = mc.simulate(ts_length, init=inits, return_values=return_values) - X_expected = np.array( - [np.arange(init, init+ts_length)%mc.n for init in inits] - ) - if return_values: - X_expected = self.state_values[X_expected] + init_idx = [1, 2] + inits = [init_idx, self.state_values[init_idx]] + path = np.array( + [np.arange(i, i+ts_length)%mc.n for i in init_idx] + ) + paths = [path, self.state_values[path]] + for method, init, X_expected in zip(methods, inits, paths): + X = getattr(mc, method)(ts_length, init) assert_array_equal(X, X_expected) From c0397df73c45bda54e93781c2d2d85a3a9b04604 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 11 Apr 2016 11:23:51 +0900 Subject: [PATCH 19/23] MarkovChain: Add test for simulate_values with init_value=None Test should fail for this commit --- quantecon/markov/tests/test_core.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 39e1c49d2..a1ed5c597 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -477,6 +477,15 @@ def test_simulate(self): X = getattr(mc, method)(ts_length, init) assert_array_equal(X, X_expected) + inits = [None, None] + seed = 1234 # init will be 2 + init_idx = 2 + path = np.arange(init_idx, init_idx+ts_length)%mc.n + paths = [path, self.state_values[path]] + for method, init, X_expected in zip(methods, inits, paths): + X = getattr(mc, method)(ts_length, init, random_state=seed) + assert_array_equal(X, X_expected) + def test_get_index(): P = [[0.4, 0.6], [0.2, 0.8]] From d761cfb968323380d8944fe206a75172d20f816b Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 11 Apr 2016 11:27:26 +0900 Subject: [PATCH 20/23] MarkovChain: Fix bug in simulate_values --- quantecon/markov/core.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index cd035e2c5..c0ee14011 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -557,7 +557,10 @@ def simulate_values(self, ts_length, init_value=None, num_reps=None, the `simulate` method for more information. """ - init_idx = self.get_index(init_value) + if init_value is not None: + init_idx = self.get_index(init_value) + else: + init_idx = None X = self.simulate(ts_length, init=init_idx, num_reps=num_reps, random_state=random_state) From 78e7e6737768ff3d80206833165b9941a5d02274 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 11 Apr 2016 23:12:22 +0900 Subject: [PATCH 21/23] DiGraph: Add *_components_indices which return indices - Change *_components to return labels - Remove get_*_components --- quantecon/graph_tools.py | 80 +++++++++++++---------------- quantecon/tests/test_graph_tools.py | 36 +++++++------ 2 files changed, 56 insertions(+), 60 deletions(-) diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py index 0e966729c..8f02c0f90 100644 --- a/quantecon/graph_tools.py +++ b/quantecon/graph_tools.py @@ -12,11 +12,11 @@ from fractions import gcd -# Decorator for get_*_components methods +# Decorator for *_components properties def annotate_nodes(func): - def new_func(self, return_labels=True): - list_of_components, return_labels = func(self, return_labels) - if return_labels and (self.node_labels is not None): + def new_func(self): + list_of_components = func(self) + if self.node_labels is not None: return [self.node_labels[c] for c in list_of_components] return list_of_components return new_func @@ -39,7 +39,7 @@ class DiGraph(object): node_labels : array_like(default=None) Array_like of length n containing the labels associated with the - states, which must be homogeneous in type. If None, the labels + nodes, which must be homogeneous in type. If None, the labels default to integers 0 through n-1. Attributes @@ -53,17 +53,27 @@ class DiGraph(object): num_strongly_connected_components : int The number of the strongly connected components. - strongly_connected_components : list(ndarray(int)) + strongly_connected_components_indices : list(ndarray(int)) List of numpy arrays containing the indices of the strongly connected components. + strongly_connected_components : list(ndarray) + List of numpy arrays containing the strongly connected + components, where the nodes are annotated with their labels (if + `node_labels` is not None). + num_sink_strongly_connected_components : int The number of the sink strongly connected components. - sink_strongly_connected_components : list(ndarray(int)) + sink_strongly_connected_components_indices : list(ndarray(int)) List of numpy arrays containing the indices of the sink strongly connected components. + sink_strongly_connected_components : list(ndarray) + List of numpy arrays containing the sink strongly connected + components, where the nodes are annotated with their labels (if + `node_labels` is not None). + is_aperiodic : bool Indicate whether the digraph is aperiodic. @@ -71,10 +81,15 @@ class DiGraph(object): The period of the digraph. Defined only for a strongly connected digraph. - cyclic_components : list(ndarray(int)) + cyclic_components_indices : list(ndarray(int)) List of numpy arrays containing the indices of the cyclic components. + cyclic_components : list(ndarray) + List of numpy arrays containing the cyclic components, where the + nodes are annotated with their labels (if `node_labels` is not + None). + References ---------- .. [1] `Strongly connected component @@ -209,28 +224,30 @@ def num_sink_strongly_connected_components(self): return len(self.sink_scc_labels) @property - def strongly_connected_components(self): + def strongly_connected_components_indices(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in range(self.num_strongly_connected_components)] + @property @annotate_nodes - def get_strongly_connected_components(self, return_labels): - return self.strongly_connected_components, return_labels + def strongly_connected_components(self): + return self.strongly_connected_components_indices @property - def sink_strongly_connected_components(self): + def sink_strongly_connected_components_indices(self): if self.is_strongly_connected: return [np.arange(self.n)] else: return [np.where(self.scc_proj == k)[0] for k in self.sink_scc_labels.tolist()] + @property @annotate_nodes - def get_sink_strongly_connected_components(self, return_labels): - return self.sink_strongly_connected_components, return_labels + def sink_strongly_connected_components(self): + return self.sink_strongly_connected_components_indices def _compute_period(self): """ @@ -303,16 +320,17 @@ def is_aperiodic(self): return (self.period == 1) @property - def cyclic_components(self): + def cyclic_components_indices(self): if self.is_aperiodic: return [np.arange(self.n)] else: return [np.where(self._cyclic_components_proj == k)[0] for k in range(self.period)] + @property @annotate_nodes - def get_cyclic_components(self, return_labels): - return self.cyclic_components, return_labels + def cyclic_components(self,): + return self.cyclic_components_indices def subgraph(self, nodes): """ @@ -342,34 +360,6 @@ def subgraph(self, nodes): return DiGraph(adj_matrix, weighted=weighted, node_labels=node_labels) -_get_method_docstr = \ -""" -Return a list of numpy arrays containing the {components}. - -Parameters ----------- -return_labels : bool(optional, default=True) - Whether to annotate the returned nodes with `node_labels`. - -Returns -------- -list(ndarray) - If `return_labels=True`, and if `node_labels` is not None, - each ndarray contains the node labels, and the node indices - (integers) otherwise. - -""" - -DiGraph.get_strongly_connected_components.__doc__ = \ - _get_method_docstr.format(components='strongly connected components') - -DiGraph.get_sink_strongly_connected_components.__doc__ = \ - _get_method_docstr.format(components='sink strongly connected components') - -DiGraph.get_cyclic_components.__doc__ = \ - _get_method_docstr.format(components='cyclic components') - - def _csr_matrix_indices(S): """ Generate the indices of nonzero entries of a csr_matrix S diff --git a/quantecon/tests/test_graph_tools.py b/quantecon/tests/test_graph_tools.py index f646e5c08..9aabccf39 100644 --- a/quantecon/tests/test_graph_tools.py +++ b/quantecon/tests/test_graph_tools.py @@ -236,16 +236,19 @@ def test_node_labels_connected_components(): sccs = [[0], [1], [2]] sink_sccs = [[0], [2]] - methods = ['get_strongly_connected_components', - 'get_sink_strongly_connected_components'] - for method, components_ind in zip(methods, [sccs, sink_sccs]): - for return_labels in [True, False]: - if return_labels: - components = [node_labels[i] for i in components_ind] - else: + properties = ['strongly_connected_components', + 'sink_strongly_connected_components'] + suffix = '_indices' + for prop0, components_ind in zip(properties, [sccs, sink_sccs]): + for return_indices in [True, False]: + if return_indices: components = components_ind + prop = prop0 + suffix + else: + components = [node_labels[i] for i in components_ind] + prop = prop0 list_of_array_equal( - sorted(getattr(g, method)(return_labels), key=lambda x: x[0]), + sorted(getattr(g, prop), key=lambda x: x[0]), sorted(components, key=lambda x: x[0]) ) @@ -257,15 +260,18 @@ def test_node_labels_cyclic_components(): cyclic_components = [[0], [1]] - methods = ['get_cyclic_components'] - for method, components_ind in zip(methods, [cyclic_components]): - for return_labels in [True, False]: - if return_labels: - components = [node_labels[i] for i in components_ind] - else: + properties = ['cyclic_components'] + suffix = '_indices' + for prop0, components_ind in zip(properties, [cyclic_components]): + for return_indices in [True, False]: + if return_indices: components = components_ind + prop = prop0 + suffix + else: + components = [node_labels[i] for i in components_ind] + prop = prop0 list_of_array_equal( - sorted(getattr(g, method)(return_labels), key=lambda x: x[0]), + sorted(getattr(g, prop), key=lambda x: x[0]), sorted(components, key=lambda x: x[0]) ) From f0c681fd37e70b3beab2a2962b602544850c70f8 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 11 Apr 2016 23:28:07 +0900 Subject: [PATCH 22/23] MarkovChain: Add *_classes_indices which return indices - Change *_classes to return values - Remove get_*_classes --- quantecon/markov/core.py | 84 ++++++++++++----------------- quantecon/markov/tests/test_core.py | 54 ++++++++----------- 2 files changed, 58 insertions(+), 80 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index c0ee14011..f69a11038 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -128,14 +128,26 @@ class MarkovChain(object): num_communication_classes : int The number of the communication classes. - communication_classes : list(ndarray(int)) - List of numpy arrays containing the communication classes. + communication_classes_indices : list(ndarray(int)) + List of numpy arrays containing the indices of the communication + classes. + + communication_classes : list(ndarray) + List of numpy arrays containing the communication classes, where + the states are annotated with their values (if `state_values` is + not None). num_recurrent_classes : int The number of the recurrent classes. - recurrent_classes : list(ndarray(int)) - List of numpy arrays containing the recurrent classes. + recurrent_classes_indices : list(ndarray(int)) + List of numpy arrays containing the indices of the recurrent + classes. + + recurrent_classes : list(ndarray) + List of numpy arrays containing the recurrent classes, where the + states are annotated with their values (if `state_values` is not + None). is_aperiodic : bool Indicate whether the Markov chain is aperiodic. @@ -143,9 +155,14 @@ class MarkovChain(object): period : int The period of the Markov chain. - cyclic_classes : list(ndarray(int)) - List of numpy arrays containing the cyclic classes. Defined only - when the Markov chain is irreducible. + cyclic_classes_indices : list(ndarray(int)) + List of numpy arrays containing the indices of the cyclic + classes. Defined only when the Markov chain is irreducible. + + cyclic_classes : list(ndarray) + List of numpy arrays containing the cyclic classes, where the + states are annotated with their values (if `state_values` is not + None). Defined only when the Markov chain is irreducible. Notes ----- @@ -315,28 +332,26 @@ def is_irreducible(self): def num_communication_classes(self): return self.digraph.num_strongly_connected_components + @property + def communication_classes_indices(self): + return self.digraph.strongly_connected_components_indices + @property def communication_classes(self): return self.digraph.strongly_connected_components - def get_communication_classes(self, return_values=True): - return self.digraph.get_strongly_connected_components( - return_labels=return_values - ) - @property def num_recurrent_classes(self): return self.digraph.num_sink_strongly_connected_components + @property + def recurrent_classes_indices(self): + return self.digraph.sink_strongly_connected_components_indices + @property def recurrent_classes(self): return self.digraph.sink_strongly_connected_components - def get_recurrent_classes(self, return_values=True): - return self.digraph.get_sink_strongly_connected_components( - return_labels=return_values - ) - @property def is_aperiodic(self): if self.is_irreducible: @@ -368,15 +383,14 @@ def cyclic_classes(self): else: return self.digraph.cyclic_components - def get_cyclic_classes(self, return_values=True): + @property + def cyclic_classes_indices(self): if not self.is_irreducible: raise NotImplementedError( 'Not defined for a reducible Markov chain' ) else: - return self.digraph.get_cyclic_components( - return_labels=return_values - ) + return self.digraph.cyclic_components_indices def _compute_stationary(self): """ @@ -652,34 +666,6 @@ def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states, out[i, t+1] = indices[indptr[out[i, t]]+k] -_get_method_docstr = \ -""" -Return a list of numpy arrays containing the {classes}. - -Parameters ----------- -return_values : bool(optional, default=True) - Whether to annotate the returned states with `state_values`. - -Returns -------- -list(ndarray) - If `return_values=True`, and if `state_values` is not None, - each ndarray contains the state values, and the state indices - (integers) otherwise. - -""" - -MarkovChain.get_communication_classes.__doc__ = \ - _get_method_docstr.format(classes='communication classes') - -MarkovChain.get_recurrent_classes.__doc__ = \ - _get_method_docstr.format(classes='recurrent classes') - -MarkovChain.get_cyclic_classes.__doc__ = \ - _get_method_docstr.format(classes='cyclic classes') - - def mc_compute_stationary(P): """ Computes stationary distributions of P, one for each recurrent diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index a1ed5c597..5dd0673a1 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -406,51 +406,43 @@ def test_com_rec_classes(self): mc = mc_dict['mc'] coms = mc_dict['coms'] recs = mc_dict['recs'] - methods = ['get_communication_classes', - 'get_recurrent_classes'] - for method, classes_ind in zip(methods, [coms, recs]): - for return_values in [True, False]: - if return_values: - classes = [self.state_values[i] for i in classes_ind] - key = lambda x: x[0, 0] - else: + properties = ['communication_classes', + 'recurrent_classes'] + suffix = '_indices' + for prop0, classes_ind in zip(properties, [coms, recs]): + for return_indices in [True, False]: + if return_indices: classes = classes_ind + prop = prop0 + suffix key = lambda x: x[0] + else: + classes = [self.state_values[i] for i in classes_ind] + prop = prop0 + key = lambda x: x[0, 0] list_of_array_equal( - sorted(getattr(mc, method)(return_values), key=key), + sorted(getattr(mc, prop), key=key), sorted(classes, key=key) ) - # Check that the default of return_values is True - classes = [self.state_values[i] for i in classes_ind] - key = lambda x: x[0, 0] - list_of_array_equal( - sorted(getattr(mc, method)(), key=key), - sorted(classes, key=key) - ) def test_cyc_classes(self): mc = self.mc_periodic_dict['mc'] cycs = self.mc_periodic_dict['cycs'] - methods = ['get_cyclic_classes'] - for method, classes_ind in zip(methods, [cycs]): - for return_values in [True, False]: - if return_values: - classes = [self.state_values[i] for i in classes_ind] - key = lambda x: x[0, 0] - else: + properties = ['cyclic_classes'] + suffix = '_indices' + for prop0, classes_ind in zip(properties, [cycs]): + for return_indices in [True, False]: + if return_indices: classes = classes_ind + prop = prop0 + suffix key = lambda x: x[0] + else: + classes = [self.state_values[i] for i in classes_ind] + prop = prop0 + key = lambda x: x[0, 0] list_of_array_equal( - sorted(getattr(mc, method)(return_values), key=key), + sorted(getattr(mc, prop), key=key), sorted(classes, key=key) ) - # Check that the default of return_values is True - classes = [self.state_values[i] for i in classes_ind] - key = lambda x: x[0, 0] - list_of_array_equal( - sorted(getattr(mc, method)(), key=key), - sorted(classes, key=key) - ) def test_simulate(self): # Deterministic mc From 05e8ad89d3adfb67618405a5df5bad5d1f90f296 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Mon, 11 Apr 2016 23:37:06 +0900 Subject: [PATCH 23/23] MarkovChain: Change simulate_values to simulate - Add simulate_indices which return indices --- quantecon/markov/core.py | 23 ++++++++++++----------- quantecon/markov/tests/test_core.py | 2 +- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py index f69a11038..1e718cbf8 100644 --- a/quantecon/markov/core.py +++ b/quantecon/markov/core.py @@ -447,9 +447,11 @@ def cdfs1d(self): self._cdfs1d = cdfs1d return self._cdfs1d - def simulate(self, ts_length, init=None, num_reps=None, random_state=None): + def simulate_indices(self, ts_length, init=None, num_reps=None, + random_state=None): """ - Simulate time series of state transitions. + Simulate time series of state transitions, where state indices + are returned. Parameters ---------- @@ -539,18 +541,17 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None): else: return X - def simulate_values(self, ts_length, init_value=None, num_reps=None, - random_state=None): + def simulate(self, ts_length, init=None, num_reps=None, random_state=None): """ - Simulate time series of state transitions, but return the values - in `state_values` instead of the state indices. + Simulate time series of state transitions, where the states are + annotated with their values (if `state_values` is not None). Parameters ---------- ts_length : scalar(int) Length of each simulation. - init_value : scalar or array_like, optional(default=None) + init : scalar or array_like, optional(default=None) Initial state values(s). If None, the initial state is randomly drawn. @@ -571,12 +572,12 @@ def simulate_values(self, ts_length, init_value=None, num_reps=None, the `simulate` method for more information. """ - if init_value is not None: - init_idx = self.get_index(init_value) + if init is not None: + init_idx = self.get_index(init) else: init_idx = None - X = self.simulate(ts_length, init=init_idx, num_reps=num_reps, - random_state=random_state) + X = self.simulate_indices(ts_length, init=init_idx, num_reps=num_reps, + random_state=random_state) # Annotate states if self.state_values is not None: diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py index 5dd0673a1..5f84a92c7 100644 --- a/quantecon/markov/tests/test_core.py +++ b/quantecon/markov/tests/test_core.py @@ -449,7 +449,7 @@ def test_simulate(self): mc = self.mc_periodic_dict['mc'] ts_length = 6 - methods = ['simulate', 'simulate_values'] + methods = ['simulate_indices', 'simulate'] init_idx = 0 inits = [init_idx, self.state_values[init_idx]]