Repository: pymc-devs/pymc
Pull request: Fixed leaf_nodes computation
This fixes #3643. I have to run now but I'll add a test, some release notes and refactor a bit of old dangling code by tomorrow. After that, it should be ok to merge.

Problem statement:
Shared vars dirty values in Deterministic RVs
sample_posterior_predictive is susceptible to "bad" ordering of var_names.
The problem arises only when using theano.shared vars for observed RVs. 
It seems that when sample_posterior_predictive is given var_names in an order in which a Deterministic RV is placed before the RV(s) it is based on, the samples will contain wrong values for that deterministic variable.
Namely - the values for the deterministic rv in the ppc trace will be values based on the "original" theano.shared values.

One possible fix is reordering var_names inside sample_posterior_predictive to be consistent with
```python
free_RVs + observed_RVs + deterministics
```
As far as I tested this ordering works, but this might be just a workaround hiding a deeper problem. Not sure.

If the reordering solution is the right way to go, better to topologically order the rvs, though I didn't find a way in the API to retrieve their dependencies.

See reproduction below.

```python
np.random.seed(123)
        obs_a = theano.shared(np.array([10., 20., 30.]))

        with pm.Model() as m:
            pm.Normal('mu', 3, 5)
            a = pm.Normal('a', 20, 10, observed=obs_a)
            pm.Deterministic('b', a * 2)

            trace = pm.sample(10)

        # var_names = ['a', 'b'] # GOOD ORDERING
        var_names = ['b', 'a'] # BAD ORDERING
        ppc = pm.sample_posterior_predictive(trace, model=m, var_names=var_names)
        sa = ppc['a'][0, 0]
        sb = ppc['b'][0, 0]
        if sa*2 != sb:
            print(f'a = {sa}')
            print(f'b = {sb}')
            raise ValueError('wrong')
        
        print('OK')
```


* PyMC3 Version: 3.7
* Theano Version: 1.0.4
* Python Version: 3.7.3
* Operating system: macOS High Sierra
* How did you install PyMC3: pip

Hints from issue discussion:
Thanks for reporting @uri-datorama. This highlights a problem in `draw_values`. In particular, in how the output from `pymc3.model.get_named_nodes_and_relations` is handled. I'll work on a fix later today.

Commit message: Fixed float32 error

Patch to review:
diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index fb2b12fa18..f8539d0646 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -28,6 +28,8 @@
 - Fixed a defect in `OrderedLogistic.__init__` that unnecessarily increased the dimensionality of the underlying `p`. Related to issue issue [#3535](https://github.com/pymc-devs/pymc3/issues/3535) but was not the true cause of it.
 - Wrapped `DensityDist.rand` with `generate_samples` to make it aware of the distribution's shape. Added control flow attributes to still be able to behave as in earlier versions, and to control how to interpret the `size` parameter in the `random` callable signature. Fixes [3553](https://github.com/pymc-devs/pymc3/issues/3553)
 - Added `theano.gof.graph.Constant` to type checks done in `_draw_value` (fixes issue [3595](https://github.com/pymc-devs/pymc3/issues/3595))
+- Refactored `pymc3.model.get_named_nodes_and_relations` to use the ancestors and descendents, in a way that is consistent with `theano`'s naming convention.
+- Changed the way in which `pymc3.model.get_named_nodes_and_relations` computes nodes without ancestors to make it robust to changes in var_name orderings (issue [#3643](https://github.com/pymc-devs/pymc3/issues/3643))
 
 
 ## PyMC3 3.7 (May 29 2019)
diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index e1188fce7a..5f2e0a375e 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -565,30 +565,30 @@ def draw_values(params, point=None, size=None):
         # specified in the point. Need to find the node-inputs, their
         # parents and children to replace them.
         leaf_nodes = {}
-        named_nodes_parents = {}
-        named_nodes_children = {}
+        named_nodes_descendents = {}
+        named_nodes_ancestors = {}
         for _, param in symbolic_params:
             if hasattr(param, 'name'):
                 # Get the named nodes under the `param` node
-                nn, nnp, nnc = get_named_nodes_and_relations(param)
+                nn, nnd, nna = get_named_nodes_and_relations(param)
                 leaf_nodes.update(nn)
                 # Update the discovered parental relationships
-                for k in nnp.keys():
-                    if k not in named_nodes_parents.keys():
-                        named_nodes_parents[k] = nnp[k]
+                for k in nnd.keys():
+                    if k not in named_nodes_descendents.keys():
+                        named_nodes_descendents[k] = nnd[k]
                     else:
-                        named_nodes_parents[k].update(nnp[k])
+                        named_nodes_descendents[k].update(nnd[k])
                 # Update the discovered child relationships
-                for k in nnc.keys():
-                    if k not in named_nodes_children.keys():
-                        named_nodes_children[k] = nnc[k]
+                for k in nna.keys():
+                    if k not in named_nodes_ancestors.keys():
+                        named_nodes_ancestors[k] = nna[k]
                     else:
-                        named_nodes_children[k].update(nnc[k])
+                        named_nodes_ancestors[k].update(nna[k])
 
         # Init givens and the stack of nodes to try to `_draw_value` from
         givens = {p.name: (p, v) for (p, size), v in drawn.items()
                   if getattr(p, 'name', None) is not None}
-        stack = list(leaf_nodes.values())  # A queue would be more appropriate
+        stack = list(leaf_nodes.values())
         while stack:
             next_ = stack.pop(0)
             if (next_, size) in drawn:
@@ -609,7 +609,7 @@ def draw_values(params, point=None, size=None):
                 # of TensorConstants or SharedVariables, we must add them
                 # to the stack or risk evaluating deterministics with the
                 # wrong values (issue #3354)
-                stack.extend([node for node in named_nodes_parents[next_]
+                stack.extend([node for node in named_nodes_descendents[next_]
                               if isinstance(node, (ObservedRV,
                                                    MultiObservedRV))
                               and (node, size) not in drawn])
@@ -618,7 +618,7 @@ def draw_values(params, point=None, size=None):
                 # If the node does not have a givens value, try to draw it.
                 # The named node's children givens values must also be taken
                 # into account.
-                children = named_nodes_children[next_]
+                children = named_nodes_ancestors[next_]
                 temp_givens = [givens[k] for k in givens if k in children]
                 try:
                     # This may fail for autotransformed RVs, which don't
@@ -633,7 +633,7 @@ def draw_values(params, point=None, size=None):
                     # The node failed, so we must add the node's parents to
                     # the stack of nodes to try to draw from. We exclude the
                     # nodes in the `params` list.
-                    stack.extend([node for node in named_nodes_parents[next_]
+                    stack.extend([node for node in named_nodes_descendents[next_]
                                   if node is not None and
                                   (node, size) not in drawn])
 
@@ -657,8 +657,8 @@ def draw_values(params, point=None, size=None):
                         # This may set values for certain nodes in the drawn
                         # dictionary, but they don't get added to the givens
                         # dictionary. Here, we try to fix that.
-                        if param in named_nodes_children:
-                            for node in named_nodes_children[param]:
+                        if param in named_nodes_ancestors:
+                            for node in named_nodes_ancestors[param]:
                                 if (
                                     node.name not in givens and
                                     (node, size) in drawn
diff --git a/pymc3/model.py b/pymc3/model.py
index 00423e5514..2c292d52e0 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -91,6 +91,7 @@ def incorporate_methods(source, destination, methods, default=None,
         else:
             setattr(destination, method, None)
 
+
 def get_named_nodes_and_relations(graph):
     """Get the named nodes in a theano graph (i.e., nodes whose name
     attribute is not None) along with their relationships (i.e., the
@@ -102,64 +103,70 @@ def get_named_nodes_and_relations(graph):
     graph - a theano node
 
     Returns:
-    leaf_nodes: A dictionary of name:node pairs, of the named nodes that
-        are also leafs of the graph
-    node_parents: A dictionary of node:set([parents]) pairs. Each key is
+    leafs: A dictionary of name:node pairs, of the named nodes that
+        have no named ancestors in the provided theano graph.
+    descendents: A dictionary of node:set([parents]) pairs. Each key is
         a theano named node, and the corresponding value is the set of
-        theano named nodes that are parents of the node. These parental
-        relations skip unnamed intermediate nodes.
-    node_children: A dictionary of node:set([children]) pairs. Each key
+        theano named nodes that are direct descendents of the node in the
+        supplied ``graph``. These relations skip unnamed intermediate nodes.
+    ancestors: A dictionary of node:set([ancestors]) pairs. Each key
         is a theano named node, and the corresponding value is the set
-        of theano named nodes that are children of the node. These child
-        relations skip unnamed intermediate nodes.
+        of theano named nodes that are direct ancestors in the of the node in
+        the supplied ``graph``. These relations skip unnamed intermediate
+        nodes.
 
     """
     if graph.name is not None:
-        node_parents = {graph: set()}
-        node_children = {graph: set()}
+        ancestors = {graph: set()}
+        descendents = {graph: set()}
     else:
-        node_parents = {}
-        node_children = {}
-    return _get_named_nodes_and_relations(graph, None, {}, node_parents, node_children)
-
-def _get_named_nodes_and_relations(graph, parent, leaf_nodes,
-                                        node_parents, node_children):
+        ancestors = {}
+        descendents = {}
+    descendents, ancestors = _get_named_nodes_and_relations(
+        graph, None, ancestors, descendents
+    )
+    leafs = {
+        node.name: node for node, ancestor in ancestors.items()
+        if len(ancestor) == 0
+    }
+    return leafs, descendents, ancestors
+
+
+def _get_named_nodes_and_relations(graph, descendent, descendents, ancestors):
     if getattr(graph, 'owner', None) is None:  # Leaf node
         if graph.name is not None:  # Named leaf node
-            leaf_nodes.update({graph.name: graph})
-            if parent is not None:  # Is None for the root node
+            if descendent is not None:  # Is None for the first node
                 try:
-                    node_parents[graph].add(parent)
+                    descendents[graph].add(descendent)
                 except KeyError:
-                    node_parents[graph] = {parent}
-                node_children[parent].add(graph)
+                    descendents[graph] = {descendent}
+                ancestors[descendent].add(graph)
             else:
-                node_parents[graph] = set()
+                descendents[graph] = set()
             # Flag that the leaf node has no children
-            node_children[graph] = set()
+            ancestors[graph] = set()
     else:  # Intermediate node
         if graph.name is not None:  # Intermediate named node
-            if parent is not None:  # Is only None for the root node
+            if descendent is not None:  # Is only None for the root node
                 try:
-                    node_parents[graph].add(parent)
+                    descendents[graph].add(descendent)
                 except KeyError:
-                    node_parents[graph] = {parent}
-                node_children[parent].add(graph)
+                    descendents[graph] = {descendent}
+                ancestors[descendent].add(graph)
             else:
-                node_parents[graph] = set()
-            # The current node will be set as the parent of the next
+                descendents[graph] = set()
+            # The current node will be set as the descendent of the next
             # nodes only if it is a named node
-            parent = graph
+            descendent = graph
             # Init the nodes children to an empty set
-            node_children[graph] = set()
+            ancestors[graph] = set()
         for i in graph.owner.inputs:
-            temp_nodes, temp_inter, temp_tree = \
-                _get_named_nodes_and_relations(i, parent, leaf_nodes,
-                                               node_parents, node_children)
-            leaf_nodes.update(temp_nodes)
-            node_parents.update(temp_inter)
-            node_children.update(temp_tree)
-    return leaf_nodes, node_parents, node_children
+            temp_desc, temp_ances = _get_named_nodes_and_relations(
+                i, descendent, descendents, ancestors
+            )
+            descendents.update(temp_desc)
+            ancestors.update(temp_ances)
+    return descendents, ancestors
 
 
 class Context:
diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py
index db431ade4f..c75a1592be 100644
--- a/pymc3/tests/test_sampling.py
+++ b/pymc3/tests/test_sampling.py
@@ -415,6 +415,27 @@ def test_deterministic_of_observed(self):
             rtol = 1e-5 if theano.config.floatX == "float64" else 1e-3
             assert np.allclose(ppc["in_1"] + ppc["in_2"], ppc["out"], rtol=rtol)
 
+    def test_var_name_order_invariance(self):
+        obs_a = theano.shared(pm.theanof.floatX(np.array([10., 20., 30.])))
+        with pm.Model() as m:
+            pm.Normal('mu', 3, 5)
+            a = pm.Normal('a', 20, 10, observed=obs_a)
+            pm.Deterministic('b', a * 2)
+            trace = pm.sample(10)
+
+        np.random.seed(123)
+        var_names = ['b', 'a']
+        ppc1 = pm.sample_posterior_predictive(
+            trace, model=m, var_names=var_names
+        )
+        np.random.seed(123)
+        var_names = ['a', 'b']
+        ppc2 = pm.sample_posterior_predictive(
+            trace, model=m, var_names=var_names
+        )
+        assert np.all(ppc1["a"] == ppc2["a"])
+        assert np.all(ppc1["b"] == ppc2["b"])
+        assert np.allclose(ppc1["b"], (2 * ppc1["a"]))
 
     def test_deterministic_of_observed_modified_interface(self):
         meas_in_1 = pm.theanof.floatX(2 + 4 * np.random.randn(100))


Code review comments:
diff_hunk: @@ -102,64 +103,70 @@ def get_named_nodes_and_relations(graph):
     graph - a theano node
 
     Returns:
-    leaf_nodes: A dictionary of name:node pairs, of the named nodes that
-        are also leafs of the graph
-    node_parents: A dictionary of node:set([parents]) pairs. Each key is
+    leafs: A dictionary of name:node pairs, of the named nodes that
+        have no named ancestors in the provided theano graph.
+    descendents: A dictionary of node:set([parents]) pairs. Each key is
         a theano named node, and the corresponding value is the set of
-        theano named nodes that are parents of the node. These parental
-        relations skip unnamed intermediate nodes.
-    node_children: A dictionary of node:set([children]) pairs. Each key
+        theano named nodes that are direct descendents of the node in the
path: pymc3/model.py
text: @user1:
Minor complaint: we say first that these are "direct descendents," but then we say "we skip intermediate nodes," so they *are not* direct descendents.  I think I know what you mean -- they are descendents with no intervening *named* nodes, but someone might be confused?

diff_hunk: @@ -102,64 +103,70 @@ def get_named_nodes_and_relations(graph):
     graph - a theano node
 
     Returns:
-    leaf_nodes: A dictionary of name:node pairs, of the named nodes that
-        are also leafs of the graph
-    node_parents: A dictionary of node:set([parents]) pairs. Each key is
+    leafs: A dictionary of name:node pairs, of the named nodes that
path: pymc3/model.py
text: @user1:
It's not terribly important but -- English is crazy -- the plural of "leaf" is "leaves" not "leafs."  We could leave (no pun intended!) it, but it's just a string replace, so we might as well fix it.

diff_hunk: @@ -565,30 +565,30 @@ def draw_values(params, point=None, size=None):
         # specified in the point. Need to find the node-inputs, their
         # parents and children to replace them.
         leaf_nodes = {}
-        named_nodes_parents = {}
-        named_nodes_children = {}
+        named_nodes_descendents = {}
path: pymc3/distributions/distribution.py
text: @user1:
While you are changing this, what about pulling this out and making it a free-standing function?  This block of code is awfully long, and I think the reader/maintainer doesn't need to know the details wouldn't it be ok to have

`named_nodes_ancestors, named_nodes_descendents = build_node_tree(symbolic_params)`