Repository: pymc-devs/pymc

Issue:
Posterior predictive doesn't resample intermediate Deterministics of intervened variables 
### Describe the issue:

When sampling with do operator, Deterministics are injected from the trace and do operation gives no effect

### Reproduceable code example:

```python
with pm.Model() as model:
    a = pm.Normal("a")
    b = pm.Deterministic("b", a + 10)
    c = pm.Deterministic("c", b * 3)
    d = pm.Normal("d", c, observed=0)
with model:
    trace = pm.sample()
with pm.do(model, {"a": -100}):
    trace1 = pm.sample_posterior_predictive(trace, var_names=["c"])
# should raise an error since a is -100
np.testing.assert_allclose(trace1.posterior_predictive["c"], trace.posterior["c"])
with pm.do(model, {"a": -100}):
    trace1 = pm.sample_posterior_predictive(trace, var_names=["c", "b"])
# raises an error since a is -100
np.testing.assert_allclose(trace1.posterior_predictive["c"], trace.posterior["c"])
```


### Error message:

```shell
...
```


### PyMC version information:

5.9.1

### Context for the issue:

_No response_

Discussion:
ricardoV94: @drbenvincent I am a bit conflicted. I do want to test that as a user, but not as part of the codebase as it's covering something that's orthogonal to it: #6876 

This PR is not a solution to that problem, but sidesteps it in this specific case. Even if we didn't have that problem we would still want this PR so you can use `set_data` if the users wished to.
drbenvincent: Sure - totally fine for you to ignore my suggestion. There are a bunch of PR's which are all closely related and I don't have enough knowledge to know which exact PR would be relevant for a user-facing test like this. Feel free to ignore. Can possibly add a test along these lines to one of the other PR's potentially.

Patch:
diff --git a/pymc/data.py b/pymc/data.py
index c21ac3001..4b3538a34 100644
--- a/pymc/data.py
+++ b/pymc/data.py
@@ -30,13 +30,12 @@ from pytensor.compile.builders import OpFromGraph
 from pytensor.compile.sharedvalue import SharedVariable
 from pytensor.graph.basic import Variable
 from pytensor.raise_op import Assert
-from pytensor.scalar import Cast
-from pytensor.tensor.elemwise import Elemwise
 from pytensor.tensor.random.basic import IntegersRV
 from pytensor.tensor.variable import TensorConstant, TensorVariable
 
 import pymc as pm
 
+from pymc.logprob.utils import rvs_in_graph
 from pymc.pytensorf import convert_data
 from pymc.vartypes import isgenerator
 
@@ -111,13 +110,7 @@ def is_valid_observed(v) -> bool:
         return True
 
     return (
-        # The only PyTensor operation we allow on observed data is type casting
-        # Although we could allow for any graph that does not depend on other RVs
-        (
-            isinstance(v.owner.op, Elemwise)
-            and isinstance(v.owner.op.scalar_op, Cast)
-            and is_valid_observed(v.owner.inputs[0])
-        )
+        not rvs_in_graph(v)
         # Or Minibatch
         or (
             isinstance(v.owner.op, MinibatchOp)
@@ -148,7 +141,7 @@ def Minibatch(variable: TensorVariable, *variables: TensorVariable, batch_size:
     for i, v in enumerate(tensors):
         if not is_valid_observed(v):
             raise ValueError(
-                f"{i}: {v} is not valid for Minibatch, only constants or constants.astype(dtype) are allowed"
+                f"{i}: {v} is not valid for Minibatch, only non-random variables are allowed"
             )
 
     upper = tensors[0].shape[0]
diff --git a/pymc/model/transform/conditioning.py b/pymc/model/transform/conditioning.py
index 7bb8abbf2..3b506d142 100644
--- a/pymc/model/transform/conditioning.py
+++ b/pymc/model/transform/conditioning.py
@@ -16,7 +16,9 @@ import warnings
 from collections.abc import Mapping, Sequence
 from typing import Any, Union
 
-from pytensor.graph import ancestors
+import pytensor
+
+from pytensor.graph import Constant, ancestors
 from pytensor.tensor import TensorVariable
 
 from pymc.logprob.transforms import Transform
@@ -126,7 +128,9 @@ def observe(
 def do(
     model: Model,
     vars_to_interventions: Mapping[Union["str", TensorVariable], Any],
-    prune_vars=False,
+    *,
+    make_interventions_shared: bool = True,
+    prune_vars: bool = False,
 ) -> Model:
     """Replace model variables by intervention variables.
 
@@ -140,6 +144,8 @@ def do(
         Dictionary that maps model variables (or names) to intervention expressions.
         Intervention expressions must have a shape and data type that is compatible
         with the original model variable.
+    make_interventions_shared: bool, defaults to True,
+        Whether to make constant interventions shared variables.
     prune_vars: bool, defaults to False
         Whether to prune model variables that are not connected to any observed variables,
         after the interventions.
@@ -170,11 +176,14 @@ def do(
 
     """
     do_mapping = {}
-    for var, obs in vars_to_interventions.items():
+    for var, intervention in vars_to_interventions.items():
         if isinstance(var, str):
             var = model[var]
         try:
-            do_mapping[var] = var.type.filter_variable(obs)
+            intervention = var.type.filter_variable(intervention)
+            if make_interventions_shared and isinstance(intervention, Constant):
+                intervention = pytensor.shared(intervention.data, name=var.name)
+            do_mapping[var] = intervention
         except TypeError as err:
             raise TypeError(
                 "Incompatible replacement type. Make sure the shape and datatype of the interventions match the original variables"
diff --git a/pymc/pytensorf.py b/pymc/pytensorf.py
index cf5700b95..82eca936b 100644
--- a/pymc/pytensorf.py
+++ b/pymc/pytensorf.py
@@ -164,6 +164,12 @@ def extract_obs_data(x: TensorVariable) -> np.ndarray:
             mask[mask_idx] = 1
             return np.ma.MaskedArray(array_data, mask)
 
+    from pymc.logprob.utils import rvs_in_graph
+
+    if not inputvars(x) and not rvs_in_graph(x):
+        cheap_eval_mode = Mode(linker="py", optimizer=None)
+        return x.eval(mode=cheap_eval_mode)
+
     raise TypeError(f"Data cannot be extracted from {x}")
 
 
diff --git a/pymc/sampling/jax.py b/pymc/sampling/jax.py
index 193554380..6c0119282 100644
--- a/pymc/sampling/jax.py
+++ b/pymc/sampling/jax.py
@@ -30,7 +30,7 @@ import pytensor.tensor as pt
 from arviz.data.base import make_attrs
 from jax.lax import scan
 from numpy.typing import ArrayLike
-from pytensor.compile import SharedVariable, Supervisor, mode
+from pytensor.compile import SharedVariable, mode
 from pytensor.graph.basic import graph_inputs
 from pytensor.graph.fg import FunctionGraph
 from pytensor.graph.replace import clone_replace
@@ -127,18 +127,6 @@ def get_jaxified_graph(
     graph = _replace_shared_variables(outputs) if outputs is not None else None
 
     fgraph = FunctionGraph(inputs=inputs, outputs=graph, clone=True)
-    # We need to add a Supervisor to the fgraph to be able to run the
-    # JAX sequential optimizer without warnings. We made sure there
-    # are no mutable input variables, so we only need to check for
-    # "destroyers". This should be automatically handled by PyTensor
-    # once https://github.com/aesara-devs/aesara/issues/637 is fixed.
-    fgraph.attach_feature(
-        Supervisor(
-            input
-            for input in fgraph.inputs
-            if not (hasattr(fgraph, "destroyers") and fgraph.has_destroyers([input]))
-        )
-    )
     mode.JAX.optimizer.rewrite(fgraph)
 
     # We now jaxify the optimized fgraph


Test patch:
diff --git a/tests/model/transform/test_conditioning.py b/tests/model/transform/test_conditioning.py
index c3dee552b..828fac737 100644
--- a/tests/model/transform/test_conditioning.py
+++ b/tests/model/transform/test_conditioning.py
@@ -16,9 +16,12 @@ import numpy as np
 import pytest
 
 from pytensor import config
+from pytensor.compile import SharedVariable
+from pytensor.graph import Constant
 
 import pymc as pm
 
+from pymc import sample_posterior_predictive, set_data
 from pymc.distributions.transforms import logodds
 from pymc.model.transform.conditioning import (
     change_value_transforms,
@@ -253,6 +256,46 @@ def test_do_self_reference():
     np.testing.assert_allclose(draw_x + 100, draw_do_x)
 
 
+def test_do_make_intervenstions_shared():
+    with pm.Model(coords={"obs": [1]}) as m:
+        x = pm.Normal("x", dims="obs")
+        y = pm.Normal("y", dims="obs")
+
+    constant_m = do(m, {x: [0.5]}, make_interventions_shared=False)
+    constant_x = constant_m["x"]
+    assert isinstance(constant_x, Constant)
+    np.testing.assert_array_equal(constant_x.data, [0.5])
+
+    shared_m = do(m, {x: [0.5]}, make_interventions_shared=True)
+    shared_x = shared_m["x"]
+    assert isinstance(shared_x, SharedVariable)
+    np.testing.assert_array_equal(shared_x.get_value(borrow=True), [0.5])
+
+    with shared_m:
+        set_data({"x": [0.6, 0.9]}, coords={"obs": [2, 3]})
+        pp_y = pm.sample_prior_predictive(draws=3).prior["y"]
+    assert pp_y.sizes == {"chain": 1, "draw": 3, "obs": 2}
+    assert pp_y.shape == (1, 3, 2)
+
+
+@pytest.mark.parametrize(
+    "make_interventions_shared",
+    [True, pytest.param(False, marks=pytest.mark.xfail(reason="#6876"))],
+)
+def test_do_sample_posterior_predictive(make_interventions_shared):
+    # Regression test for https://github.com/pymc-devs/pymc/issues/6977
+    with pm.Model() as model:
+        a = pm.Normal("a")
+        b = pm.Deterministic("b", a * 2)
+        c = pm.Normal("c", b / 2)
+
+    idata = az.from_dict({"a": [[1.0]], "b": [[2.0]], "c": [[1.0]]})
+
+    with do(model, {a: 1000}, make_interventions_shared=make_interventions_shared):
+        pp = sample_posterior_predictive(idata, var_names=["c"], predictions=True).predictions
+    assert (pp["c"] > 500).all()
+
+
 def test_change_value_transforms():
     with pm.Model() as base_m:
         p = pm.Uniform("p", 0, 1, default_transform=None)
diff --git a/tests/test_data.py b/tests/test_data.py
index 154737b63..afca1831a 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -509,11 +509,17 @@ class TestMinibatch:
         mb = pm.Minibatch(pt.as_tensor(self.data).astype(int), batch_size=20)
         assert isinstance(mb.owner.op, MinibatchOp)
 
-        with pytest.raises(ValueError, match="not valid for Minibatch"):
-            pm.Minibatch(pt.as_tensor(self.data) * 2, batch_size=20)
+        mb = pm.Minibatch(pt.as_tensor(self.data) * 2, batch_size=20)
+        assert isinstance(mb.owner.op, MinibatchOp)
+
+        for mb in pm.Minibatch(self.data, pt.as_tensor(self.data) * 2, batch_size=20):
+            assert isinstance(mb.owner.op, MinibatchOp)
 
-        with pytest.raises(ValueError, match="not valid for Minibatch"):
-            pm.Minibatch(self.data, pt.as_tensor(self.data) * 2, batch_size=20)
+    def test_not_allowed(self):
+        data = pt.random.normal(loc=self.data, scale=1)
+
+        with pytest.raises(ValueError):
+            pm.Minibatch(data, batch_size=20)
 
     def test_assert(self):
         d1, d2 = pm.Minibatch(self.data, self.data[::2], batch_size=20)
@@ -530,3 +536,21 @@ class TestMinibatch:
         [draw_mA, draw_mB] = pm.draw([mA, mB])
         assert draw_mA.shape == (10,)
         np.testing.assert_allclose(draw_mA, -draw_mB)
+
+
+def test_scaling_data_works_in_likelihood() -> None:
+    data = np.array([10, 11, 12, 13, 14, 15])
+
+    with pm.Model():
+        target = pm.Data("target", data)
+        scale = 12
+        scaled_target = target / scale
+        mu = pm.Normal("mu", mu=0, sigma=1)
+        pm.Normal("x", mu=mu, sigma=1, observed=scaled_target)
+
+        idata = pm.sample(10, chains=1, tune=100)
+
+    np.testing.assert_allclose(
+        idata.observed_data["x"].values,
+        data / scale,
+    )
diff --git a/tests/test_pytensorf.py b/tests/test_pytensorf.py
index b8c82886b..34360397a 100644
--- a/tests/test_pytensorf.py
+++ b/tests/test_pytensorf.py
@@ -195,6 +195,21 @@ class TestExtractObsData:
         assert isinstance(res, np.ndarray)
         np.testing.assert_array_equal(res, y)
 
+    def test_pytensor_operations(self):
+        x = np.array([1, 2, 3])
+        target = 1 + 3 * pt.as_tensor_variable(x)
+
+        res = extract_obs_data(target)
+        assert isinstance(res, np.ndarray)
+        np.testing.assert_array_equal(res, np.array([4, 7, 10]))
+
+    def test_pytensor_operations_raises(self):
+        x = pt.scalar("x")
+        target = 1 + 3 * x
+
+        with pytest.raises(TypeError, match="Data cannot be extracted from"):
+            extract_obs_data(target)
+
 
 @pytest.mark.parametrize("input_dtype", ["int32", "int64", "float32", "float64"])
 def test_convert_data(input_dtype):