diff --git a/doc/modules/high_level_api.rst b/doc/modules/high_level_api.rst
index bde25aa3..ca832934 100644
--- a/doc/modules/high_level_api.rst
+++ b/doc/modules/high_level_api.rst
@@ -38,6 +38,7 @@ Models
 
    ConceptBottleneckModel
    ConceptEmbeddingModel
+   ConceptMemoryReasoner
    GraphConceptBottleneckModel
    CausallyReliableConceptBottleneckModel
    BlackBox
diff --git a/examples/contributing/model.md b/examples/contributing/model.md
index be328f4f..5dbebe34 100644
--- a/examples/contributing/model.md
+++ b/examples/contributing/model.md
@@ -393,7 +393,6 @@ from torch_concepts.nn import (
 #### Special Layers
 ```python
 from torch_concepts.nn import (
-    SelectorLatentToExogenous,          # Memory-augmented selection
     WANDAGraphLearner,       # Learn concept graph structure
 )
 ```
diff --git a/examples/utilization/0_layer/7_concept_based_memory_reasoner.py b/examples/utilization/0_layer/7_concept_based_memory_reasoner.py
new file mode 100644
index 00000000..8b085a18
--- /dev/null
+++ b/examples/utilization/0_layer/7_concept_based_memory_reasoner.py
@@ -0,0 +1,115 @@
+"""
+Example: Concept Memory Reasoner with Low-Level API
+
+This example demonstrates how to build a Concept Memory Reasoner (CMR)
+using the low-level encoder and predictor layers.
+"""
+import torch
+from sklearn.metrics import accuracy_score
+from torch.nn import ModuleDict
+
+from torch_concepts import seed_everything
+from torch_concepts.data.datasets import ToyDataset
+from torch_concepts.nn import (
+    RuleMemory,
+    RuleReconstructionPredictor,
+    RuleTaskPredictor,
+    CategoricalSelector,
+    LinearEmbeddingToConcept,
+)
+
+
+def main():
+    latent_dims = 10
+    n_epochs = 500
+    n_samples = 1000
+    nb_rules = 10
+    memory_latent_size = 100
+
+    seed_everything(42)
+
+    dataset = ToyDataset(dataset='xor', seed=42, n_gen=n_samples)
+    x_train = dataset.input_data
+    concept_idx = list(dataset.graph.edge_index[0].unique().numpy())
+    task_idx = list(dataset.graph.edge_index[1].unique().numpy())
+    c_train = dataset.concepts[:, concept_idx]
+    y_train = dataset.concepts[:, task_idx]
+
+    n_features = x_train.shape[1]
+    n_concepts = c_train.shape[1]
+    n_tasks = y_train.shape[1]
+
+    latent_encoder = torch.nn.Sequential(
+        torch.nn.Linear(n_features, latent_dims),
+        torch.nn.LeakyReLU(),
+    )
+    selector_encoder = CategoricalSelector(
+        in_latent=latent_dims,
+        out_concepts=n_tasks,
+        out_exogenous=nb_rules,
+    )
+    concept_encoder = LinearEmbeddingToConcept(in_embeddings=latent_dims, out_concepts=n_concepts)
+    memory = RuleMemory(
+        n_tasks=n_tasks,
+        n_rules=nb_rules,
+        n_concepts=n_concepts,
+        latent_size=memory_latent_size,
+    )
+    task_predictor = RuleTaskPredictor(
+        in_concepts=n_concepts,
+        in_exogenous=nb_rules,
+        out_concepts=n_tasks,
+    )
+    reconstruction_predictor = RuleReconstructionPredictor(
+        in_concepts=n_concepts,
+        in_exogenous=nb_rules,
+        out_concepts=n_tasks,
+        rec_weight=0.1,
+    )
+
+    model = ModuleDict({
+        'latent_encoder': latent_encoder,
+        'selector_encoder': selector_encoder,
+        'concept_encoder': concept_encoder,
+        'memory': memory,
+        'task_predictor': task_predictor,
+        'reconstruction_predictor': reconstruction_predictor,
+    })
+
+    optimizer = torch.optim.AdamW(model.parameters(), lr=0.01)
+    concept_loss_fn = torch.nn.BCEWithLogitsLoss()
+    task_loss_fn = torch.nn.BCELoss(reduction='none')
+    model.train()
+
+    for epoch in range(n_epochs):
+        optimizer.zero_grad()
+
+        emb = latent_encoder(x_train)
+        selector = selector_encoder(latent=emb)
+        c_logits = concept_encoder(embeddings=emb)
+        c_probs = c_logits.sigmoid()
+        roles = memory()
+
+        y_pred = task_predictor(concepts=c_probs, selector=selector, roles=roles)
+        y_pred_with_rec = reconstruction_predictor(concepts=c_probs, selector=selector, roles=roles)
+
+        concept_loss = concept_loss_fn(c_logits, c_train)
+        task_loss_no_rec = task_loss_fn(y_pred, y_train)
+        task_loss_with_rec = task_loss_fn(y_pred_with_rec, y_train)
+        switched_task_loss = ((1.0 - y_train) * task_loss_no_rec + y_train * task_loss_with_rec).mean()
+        loss = concept_loss + switched_task_loss
+
+        loss.backward()
+        optimizer.step()
+
+        if epoch % 100 == 0:
+            task_accuracy = accuracy_score(y_train.cpu(), (y_pred.detach() > 0.5).cpu())
+            concept_accuracy = accuracy_score(c_train.cpu(), (c_logits.detach() > 0.0).cpu())
+            print(
+                f'Epoch {epoch}: Loss {loss.item():.2f} | '
+                f'Task Acc: {task_accuracy:.2f} | Concept Acc: {concept_accuracy:.2f}'
+            )
+
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/utilization/2.2_model/10_different_training_modes.py b/examples/utilization/2.2_model/10_different_training_modes.py
index a7146875..5dd9c64f 100644
--- a/examples/utilization/2.2_model/10_different_training_modes.py
+++ b/examples/utilization/2.2_model/10_different_training_modes.py
@@ -22,9 +22,10 @@
 
 import torch
 from torch_concepts import seed_everything
-from torch_concepts.nn import ConceptBottleneckModel, ConceptEmbeddingModel, MLP
+from torch_concepts.nn import ConceptBottleneckModel, ConceptEmbeddingModel, MLP, ConceptMemoryReasoner, CMRBlendedLoss, ConceptLoss
 from torch_concepts.nn import DeterministicInference, IndependentInference
 from torch_concepts.data import ToyDataset
+from torch_concepts.data.datasets import ToyDataset
 from torch_concepts.data.base.datamodule import ConceptDataModule
 from torch.distributions import Bernoulli
 
@@ -34,7 +35,7 @@
 
 
 def evaluate(model, datamodule, n_concepts, query):
-    """Evaluate model on test set and return concept/task accuracy."""
+    """Evaluate model on a data split and return concept/task accuracy."""
     concept_acc_fn = BinaryAccuracy()
     task_acc_fn = BinaryAccuracy()
 
@@ -48,10 +49,12 @@ def evaluate(model, datamodule, n_concepts, query):
         for batch in test_loader:
             # model.eval() automatically selects eval_inference
             out = model(input=batch['inputs']['x'], query=query)
-            c_logits = out.logits[:, :n_concepts]
-            y_logits = out.logits[:, n_concepts:]
-            c_pred = torch.sigmoid(c_logits)
-            y_pred = torch.sigmoid(y_logits)
+            predictions = out.logits if out.logits is not None else out.probs
+            c_pred = predictions[:, :n_concepts]
+            y_pred = predictions[:, n_concepts:]
+            if out.logits is not None:
+                c_pred = torch.sigmoid(c_pred)
+                y_pred = torch.sigmoid(y_pred)
 
             c_true = batch['concepts']['c'][:, :n_concepts]
             y_true = batch['concepts']['c'][:, n_concepts:]
@@ -104,7 +107,7 @@ def main():
 
     # Define variable distributions as Bernoulli
     variable_distributions = {name: Bernoulli for name in concept_names}
-    loss = torch.nn.BCEWithLogitsLoss()
+    loss = ConceptLoss(binary=torch.nn.BCEWithLogitsLoss(), binary_param="logits")
     optim = torch.optim.AdamW
     optim_kwargs = {'lr': 0.1}
 
@@ -202,6 +205,45 @@ def main():
     trainer_cem.fit(model_cem, datamodule=datamodule)
     evaluate(model_cem, datamodule, n_concepts, query)
 
+    # =========================================================================
+    # CMR WITH JOINT TRAINING
+    # =========================================================================
+    print("\n" + "=" * 60)
+    print("Example 4: CMR with Joint Training")
+    print("=" * 60)
+    print("Uses DeterministicInference for both training and evaluation")
+
+    cmr_loss = CMRBlendedLoss(task_names=['xor'])
+    optim_kwargs_cmr = {'lr': 0.01}
+
+    model_cmr = ConceptMemoryReasoner(
+        input_size=n_features,
+        annotations=annotations,
+        backbone=MLP(input_size=n_features, hidden_size=16, n_layers=1),
+        latent_size=16,
+        variable_distributions=variable_distributions,
+        task_names=['xor'],
+        n_rules=10,
+        memory_latent_size=100,
+        memory_decoder_hidden_layers=1,
+        selector_hidden_layers=1,
+        hard_roles_at_eval=True,
+        inference=DeterministicInference,
+        train_inference=DeterministicInference,
+        lightning=True,
+        loss=cmr_loss,
+        rec_weight=0,
+        optim_class=optim,
+        optim_kwargs=optim_kwargs_cmr,
+    )
+    print(f"Model type: {type(model_cmr).__name__}")
+    print(f"Eval inference: {model_cmr.eval_inference.__class__.__name__}")
+    print(f"Training inference: {model_cmr.train_inference.__class__.__name__}")
+
+    trainer_cmr = Trainer(max_epochs=100)
+    trainer_cmr.fit(model_cmr, datamodule=datamodule)
+    evaluate(model_cmr, datamodule, n_concepts, query)
+
 
 if __name__ == "__main__":
     main()
\ No newline at end of file
diff --git a/tests/nn/modules/high/models/test_cmr.py b/tests/nn/modules/high/models/test_cmr.py
new file mode 100644
index 00000000..209bcb68
--- /dev/null
+++ b/tests/nn/modules/high/models/test_cmr.py
@@ -0,0 +1,26 @@
+import torch
+from torch_concepts import Annotations
+from torch_concepts.nn import CMRBlendedLoss
+from torch_concepts.nn.modules.high.models.cmr import ConceptMemoryReasoner
+
+
+def test_cmr_routes_reconstruction_prediction_through_modeloutput_extra():
+    model = ConceptMemoryReasoner(
+        input_size=2,
+        annotations=Annotations(labels=["c1", "c2", "xor"], cardinalities=[1, 1, 1]),
+        task_names=["xor"],
+        n_rules=3,
+    )
+    target = torch.tensor([[0., 1., 1.], [1., 0., 0.]])
+    query = model.build_query(target)
+    query["tasks_with_rec"] = None
+    output = model(query=query, evidence={"input": torch.randn(2, 2)})
+
+    assert output.probs["xor"].shape == (2, 1)
+    assert output.extra["task_input"].shape == (2, 1)
+    assert output.extra["input_with_rec"].shape == (2, 1)
+    assert "tasks_with_rec" not in output.probs.annotation.label_to_index
+
+    loss = CMRBlendedLoss(task_names=["xor"])(output, model.prepare_target(target))
+    loss.backward()
+    assert torch.isfinite(loss)
diff --git a/tests/nn/modules/low/encoders/test_selector.py b/tests/nn/modules/low/encoders/test_selector.py
index 3aef7e7c..b3d59ced 100644
--- a/tests/nn/modules/low/encoders/test_selector.py
+++ b/tests/nn/modules/low/encoders/test_selector.py
@@ -3,6 +3,7 @@
 import torch
 import torch.nn as nn
 from torch_concepts.nn.modules.low.dense_layers import SelectorEmbeddingEncoder
+from torch_concepts.nn.modules.low.encoders.selector import CategoricalSelector
 
 
 class TestSelectorEmbeddingEncoder(unittest.TestCase):
@@ -119,5 +120,111 @@ def test_batch_processing(self):
             self.assertEqual(output.shape, (batch_size, 3, 4))
 
 
+class TestCategoricalSelector(unittest.TestCase):
+    """Test CategoricalSelector."""
+
+    def test_initialization(self):
+        """Test selector initialization."""
+        selector = CategoricalSelector(
+            in_latent=64,
+            out_concepts=5,
+            out_exogenous=8,
+            selector_hidden_layers=2,
+        )
+        self.assertEqual(selector.in_latent, 64)
+        self.assertEqual(selector.out_concepts, 5)
+        self.assertEqual(selector.out_exogenous, 8)
+        self.assertEqual(selector.selector_hidden_layers, 2)
+
+    def test_forward_shape(self):
+        """Test forward pass output shape."""
+        selector = CategoricalSelector(
+            in_latent=64,
+            out_concepts=4,
+            out_exogenous=6,
+        )
+        latent = torch.randn(2, 64)
+        output = selector(latent=latent)
+        self.assertEqual(output.shape, (2, 4, 6))
+
+    def test_output_is_normalized_over_exogenous_dim(self):
+        """Test output probabilities sum to 1 over exogenous dimension."""
+        selector = CategoricalSelector(
+            in_latent=32,
+            out_concepts=3,
+            out_exogenous=5,
+        )
+        latent = torch.randn(3, 32)
+        output = selector(latent=latent)
+
+        sums = output.sum(dim=-1)
+        self.assertTrue(torch.allclose(sums, torch.ones_like(sums), atol=1e-5))
+
+    def test_gradient_flow(self):
+        """Test gradient flow through selector."""
+        selector = CategoricalSelector(
+            in_latent=32,
+            out_concepts=3,
+            out_exogenous=4,
+        )
+        embeddings = torch.randn(2, 32, requires_grad=True)
+        output = selector(latent=embeddings)
+        loss = output.sum()
+        loss.backward()
+        self.assertIsNotNone(embeddings.grad)
+
+    def test_hidden_layers_configuration(self):
+        """Test configurable hidden layers in selector network."""
+        selector_zero = CategoricalSelector(
+            in_latent=32,
+            out_concepts=3,
+            out_exogenous=4,
+            selector_hidden_layers=0,
+        )
+        selector_two = CategoricalSelector(
+            in_latent=32,
+            out_concepts=3,
+            out_exogenous=4,
+            selector_hidden_layers=2,
+        )
+
+        linear_zero = sum(isinstance(layer, nn.Linear) for layer in selector_zero.selector)
+        linear_two = sum(isinstance(layer, nn.Linear) for layer in selector_two.selector)
+
+        self.assertEqual(linear_zero, 1)
+        self.assertEqual(linear_two, 3)
+
+    def test_selector_hidden_layers_validation(self):
+        """Test hidden layer argument validation."""
+        with self.assertRaises(ValueError):
+            CategoricalSelector(
+                in_latent=32,
+                out_concepts=3,
+                out_exogenous=4,
+                selector_hidden_layers=-1,
+            )
+
+    def test_selector_network(self):
+        """Test selector network structure."""
+        selector = CategoricalSelector(
+            in_latent=64,
+            out_concepts=4,
+            out_exogenous=6,
+        )
+        self.assertIsInstance(selector.selector, nn.Sequential)
+
+    def test_batch_processing(self):
+        """Test different batch sizes."""
+        selector = CategoricalSelector(
+            in_latent=32,
+            out_concepts=3,
+            out_exogenous=4,
+        )
+        for batch_size in [1, 4, 8]:
+            embeddings = torch.randn(batch_size, 32)
+            output = selector(latent=embeddings)
+            self.assertEqual(output.shape, (batch_size, 3, 4))
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/nn/modules/low/predictors/test_rules.py b/tests/nn/modules/low/predictors/test_rules.py
new file mode 100644
index 00000000..aafb4c96
--- /dev/null
+++ b/tests/nn/modules/low/predictors/test_rules.py
@@ -0,0 +1,93 @@
+"""Comprehensive tests for torch_concepts.nn.modules.low.predictors."""
+import unittest
+
+import torch
+
+from torch_concepts.nn import (
+    RuleMemory,
+    RuleReconstructionPredictor,
+    RuleTaskPredictor,
+)
+
+class TestRuleMemory(unittest.TestCase):
+    def test_initialization(self):
+        memory = RuleMemory(n_tasks=3, n_rules=5, n_concepts=10, latent_size=64, hidden_layers=2)
+        self.assertEqual(memory.shape, (3, 5, 10, 3))
+        self.assertEqual(memory.memory.weight.shape, (3, 64))
+
+    def test_forward_shape(self):
+        memory = RuleMemory(n_tasks=2, n_rules=4, n_concepts=6)
+        roles = memory()
+        self.assertEqual(roles.shape, (2, 4, 6, 3))
+        self.assertTrue(torch.all((roles >= 0) & (roles <= 1)))
+        self.assertTrue(torch.allclose(roles.sum(dim=-1), torch.ones_like(roles.sum(dim=-1))))
+
+    def test_hidden_layers_config(self):
+        memory_zero = RuleMemory(n_tasks=2, n_rules=3, n_concepts=4, hidden_layers=0)
+        memory_two = RuleMemory(n_tasks=2, n_rules=3, n_concepts=4, hidden_layers=2)
+        linear_zero = sum(isinstance(layer, torch.nn.Linear) for layer in memory_zero.decoder)
+        linear_two = sum(isinstance(layer, torch.nn.Linear) for layer in memory_two.decoder)
+        self.assertEqual(linear_zero, 1)
+        self.assertEqual(linear_two, 3)
+
+    def test_gradient_flow(self):
+        memory = RuleMemory(n_tasks=2, n_rules=3, n_concepts=4)
+        loss = memory().sum()
+        loss.backward()
+        self.assertIsNotNone(memory.memory.weight.grad)
+
+
+class TestRuleTaskPredictor(unittest.TestCase):
+    def test_forward_shape(self):
+        predictor = RuleTaskPredictor(in_concepts=6, in_exogenous=3, out_concepts=2)
+        concepts = torch.rand(4, 6)
+        selector = torch.softmax(torch.randn(4, 2, 3), dim=-1)
+        roles = torch.softmax(torch.randn(2, 3, 6, 3), dim=-1)
+        output = predictor(concepts=concepts, selector=selector, roles=roles)
+        self.assertEqual(output.shape, (4, 2))
+
+    def test_gradient_flow_detaches_concepts(self):
+        predictor = RuleTaskPredictor(in_concepts=5, in_exogenous=4, out_concepts=2)
+        concepts = torch.rand(2, 5, requires_grad=True)
+        selector = torch.softmax(torch.randn(2, 2, 4), dim=-1).requires_grad_()
+        roles = torch.softmax(torch.randn(2, 4, 5, 3), dim=-1).requires_grad_()
+        output = predictor(concepts=concepts, selector=selector, roles=roles)
+        output.sum().backward()
+        self.assertIsNone(concepts.grad)
+        self.assertIsNotNone(selector.grad)
+        self.assertIsNotNone(roles.grad)
+
+
+class TestRuleReconstructionPredictor(unittest.TestCase):
+    def test_forward_shape(self):
+        predictor = RuleReconstructionPredictor(in_concepts=6, in_exogenous=3, out_concepts=2, rec_weight=0.5)
+        concepts = torch.rand(4, 6)
+        selector = torch.softmax(torch.randn(4, 2, 3), dim=-1)
+        roles = torch.softmax(torch.randn(2, 3, 6, 3), dim=-1)
+        output = predictor(concepts=concepts, selector=selector, roles=roles)
+        self.assertEqual(output.shape, (4, 2))
+
+    def test_rec_weight_changes_output(self):
+        concepts = torch.rand(3, 4)
+        selector = torch.softmax(torch.randn(3, 2, 2), dim=-1)
+        roles = torch.softmax(torch.randn(2, 2, 4, 3), dim=-1)
+        low = RuleReconstructionPredictor(in_concepts=4, in_exogenous=2, out_concepts=2, rec_weight=0.0)
+        high = RuleReconstructionPredictor(in_concepts=4, in_exogenous=2, out_concepts=2, rec_weight=1.0)
+        out_low = low(concepts=concepts, selector=selector, roles=roles)
+        out_high = high(concepts=concepts, selector=selector, roles=roles)
+        self.assertFalse(torch.allclose(out_low, out_high))
+
+    def test_gradient_flow_detaches_concepts(self):
+        predictor = RuleReconstructionPredictor(in_concepts=5, in_exogenous=4, out_concepts=2)
+        concepts = torch.rand(2, 5, requires_grad=True)
+        selector = torch.softmax(torch.randn(2, 2, 4), dim=-1).requires_grad_()
+        roles = torch.softmax(torch.randn(2, 4, 5, 3), dim=-1).requires_grad_()
+        output = predictor(concepts=concepts, selector=selector, roles=roles)
+        output.sum().backward()
+        self.assertIsNone(concepts.grad)
+        self.assertIsNotNone(selector.grad)
+        self.assertIsNotNone(roles.grad)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/nn/modules/test_loss.py b/tests/nn/modules/test_loss.py
index ac39effb..e843e185 100644
--- a/tests/nn/modules/test_loss.py
+++ b/tests/nn/modules/test_loss.py
@@ -1301,4 +1301,4 @@ def test_single_categorical_no_padding(self):
 
 
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()
\ No newline at end of file
diff --git a/torch_concepts/data/datasets/__init__.py b/torch_concepts/data/datasets/__init__.py
index 655a0a91..87b6967d 100644
--- a/torch_concepts/data/datasets/__init__.py
+++ b/torch_concepts/data/datasets/__init__.py
@@ -1 +1,5 @@
-__all__: list[str] = []
\ No newline at end of file
+"""Built-in concept datasets."""
+
+from .toy import ToyDataset
+
+__all__ = ["ToyDataset"]
diff --git a/torch_concepts/nn/__init__.py b/torch_concepts/nn/__init__.py
index 084f8258..9dc31aac 100644
--- a/torch_concepts/nn/__init__.py
+++ b/torch_concepts/nn/__init__.py
@@ -27,11 +27,17 @@
 from .modules.low.encoders.linear import LinearEmbeddingToConcept
 from .modules.low.encoders.whitening import ConceptWhitening, WhitenedEmbeddingToConcept
 from .modules.low.encoders.cav import CAVEmbeddingToConcept
+from .modules.low.encoders.selector import CategoricalSelector
 
 # Predictors
 from .modules.low.predictors.call import CallableConceptToConcept
 from .modules.low.predictors.hypernet import HyperlinearConceptEmbeddingToConcept
 from .modules.low.predictors.linear import LinearConceptToConcept
+from .modules.low.predictors.rule import (
+    RuleMemory,
+    RuleTaskPredictor,
+    RuleReconstructionPredictor,
+)
 from .modules.low.predictors.mix import MixConceptEmbeddingToConcept
 
 # Dense layers
@@ -43,7 +49,7 @@
 
 # Loss functions
 from .modules.loss import ConceptLoss, WeightedConceptLoss, DepthWeightedConceptLoss, \
-    L1LogitRegularizer
+    L1LogitRegularizer, CMRBlendedLoss
 
 # Metrics
 from .modules.metrics import ConceptMetrics, compute_cace
@@ -57,6 +63,7 @@
 from .modules.high.models.cem import ConceptEmbeddingModel
 from .modules.high.models.graph_cbm import GraphConceptBottleneckModel
 from .modules.high.models.c2bm import CausallyReliableConceptBottleneckModel
+from .modules.high.models.cmr import ConceptMemoryReasoner
 
 # Models (mid-level)
 from .modules.mid.factors.factor import ParametricFactor
@@ -130,6 +137,9 @@
 
     # Predictor classes
     "LinearConceptToConcept",
+    "RuleMemory",
+    "RuleTaskPredictor",
+    "RuleReconstructionPredictor",
     "CallableConceptToConcept",
     "HyperlinearConceptEmbeddingToConcept",
     "MixConceptEmbeddingToConcept",
@@ -138,6 +148,8 @@
     "Dense",
     "ResidualMLP",
     "MLP",
+
+    "CategoricalSelector",
     "Sequential",
     "LinearEmbeddingEncoder",
     "SelectorEmbeddingEncoder",
@@ -146,6 +158,7 @@
     "WANDAGraphLearner",
 
     # Loss functions
+    "CMRBlendedLoss",
     "ConceptLoss",
     "WeightedConceptLoss",
     "DepthWeightedConceptLoss",
@@ -164,6 +177,7 @@
     "BlackBoxTaskOnly",
     "ConceptBottleneckModel",
     "ConceptEmbeddingModel",
+    "ConceptMemoryReasoner",
     "GraphConceptBottleneckModel",
     "CausallyReliableConceptBottleneckModel",
 
diff --git a/torch_concepts/nn/modules/high/models/cmr.py b/torch_concepts/nn/modules/high/models/cmr.py
new file mode 100644
index 00000000..419dca73
--- /dev/null
+++ b/torch_concepts/nn/modules/high/models/cmr.py
@@ -0,0 +1,110 @@
+"""Concept-based Memory Reasoner (CMR)."""
+from typing import List, Optional, Union
+import torch
+from torch import nn
+from torch.distributions import Bernoulli, OneHotCategorical
+from .....annotations import Annotations
+from .....distributions import Delta
+from ...low.encoders.linear import LinearEmbeddingToConcept
+from ...low.encoders.selector import CategoricalSelector
+from ...low.predictors.rule import RuleMemory, RuleTaskPredictor, RuleReconstructionPredictor
+from ...low.priors import LearnablePrior
+from ...mid.distributions import DEFAULT_DIST_KWARGS
+from ...mid.factors.cpd import ParametricCPD
+from ...mid.graph.bayesian_network import BayesianNetwork
+from ...mid.inference.base import BaseInference
+from ...mid.inference.torch.deterministic import DeterministicInference
+from ...mid.variable import EmbeddingVariable
+from ...outputs import ModelOutput
+from ..base.bipartite import BipartiteModel
+
+class ConceptMemoryReasoner(BipartiteModel):
+    """A neurosymbolic concept-based models that performs gradient-based rule learning and instance-wise rule selection.
+
+    The model keeps the standard bipartite concept-task structure but introduces
+    three latent objects inside the PGM:
+
+    - rule_selector: per-task categorical rule weights predicted from the
+      latent representation;
+    - rule_roles: decoded concept role probabilities from the
+      learned memory, representing propositional logic rules;
+    - tasks_with_rec: an auxiliary reconstruction output based on the roles that can be used for regularization.
+
+    Args:
+        input_size: Number of input features.
+        annotations: Dataset annotations containing binary concepts and tasks.
+        task_names: Name or names of the task labels.
+        n_rules: Number of rules stored per task.
+        memory_latent_size: Size of each learned task memory embedding.
+        memory_decoder_hidden_layers: Number of hidden decoder layers in ``RuleMemory``.
+        selector_hidden_layers: Number of hidden layers in the rule selector.
+        rec_weight: Non-negative exponent applied to each rule reconstruction probability.
+        hard_roles_at_eval: If true, use the argmax one-hot role assignment in evaluation mode.
+
+    References:
+        Debot et al. "Interpretable Concept-Based Memory Reasoning", NeurIPS 2024.
+        https://arxiv.org/abs/2407.15527
+    """
+    supported_concept_types = frozenset({"binary"})
+    param_for_discrete_var = "probs"
+    variable_distributions = {"binary": Bernoulli}
+    variable_dist_kwargs = dict(DEFAULT_DIST_KWARGS)
+    def __init__(self, input_size: int, annotations: Annotations, task_names: Union[List[str], str], n_rules: int = 10, memory_latent_size: int = 100, memory_decoder_hidden_layers: int = 1, selector_hidden_layers: int = 1, rec_weight: float = 0.1, hard_roles_at_eval: bool = True, inference: Optional[BaseInference] = DeterministicInference, inference_kwargs: Optional[dict] = None, train_inference: Optional[BaseInference] = None, train_inference_kwargs: Optional[dict] = None, lightning: bool = False, plate: Optional[bool] = None, **kwargs):
+        super().__init__(input_size=input_size, annotations=annotations, task_names=task_names, lightning=lightning, plate=plate, **kwargs)
+        if any(self.concept_annotations.concept(n).cardinality != 1 for n in self.intermediate_concept_names):
+            raise ValueError("ConceptMemoryReasoner requires binary scalar concepts.")
+        self.n_rules, self.memory_latent_size = n_rules, memory_latent_size
+        self.memory_decoder_hidden_layers, self.selector_hidden_layers = memory_decoder_hidden_layers, selector_hidden_layers
+        self.rec_weight = rec_weight
+        self.hard_roles_at_eval = hard_roles_at_eval
+        self.pgm = self._build_model()
+        self.setup_inference(inference, inference_kwargs, train_inference, train_inference_kwargs)
+
+    def default_query(self, ground_truth):
+        """Train both CMR task paths in one inference query."""
+        query = super().default_query(ground_truth)
+        query["tasks_with_rec"] = None
+        return query
+
+    def forward(self, query, evidence=None, input=None, **inference_kwargs):
+        out = super().forward(query=query, evidence=evidence, input=input, **inference_kwargs)
+        probs = out.probs
+        if probs is not None and "tasks_with_rec" in probs.annotation.label_to_index:
+            rec = probs["tasks_with_rec"]
+            keep = [name for name in self.concept_names if name in probs.annotation.label_to_index]
+            ordinary = probs[self.task_names]
+            out.params["probs"] = probs[keep]
+            extra = dict(out.extra) if out.extra else {}
+            extra["task_input"] = ordinary
+            extra["input_with_rec"] = rec
+            out.extra = extra
+        return out
+
+    def _build_model(self) -> BayesianNetwork:
+        input_var = EmbeddingVariable("input", distribution=Delta, shape=self.input_size)
+        input_cpd = ParametricCPD(input_var, parents=[], parametrization={"value": LearnablePrior(input_var.shape)})
+
+        latent_var = EmbeddingVariable("latent", distribution=Delta, size=self.latent_size)
+        latent_cpd = ParametricCPD(latent_var, parents=[input_var], parametrization={"value": self.backbone})
+
+        concepts = self.build_concept_variables(self.intermediate_concept_names, "concepts")
+        concept_cpds = ParametricCPD(concepts, parents=[latent_var], parametrization=[{"probs": nn.Sequential(LinearEmbeddingToConcept(self.latent_size, c.size), nn.Sigmoid())} for c in concepts])
+        n_concepts = sum(c.size for c in concepts)
+
+        selector = EmbeddingVariable("rule_selector", distribution=OneHotCategorical, shape=(len(self.task_names), self.n_rules))
+        selector_cpd = ParametricCPD(selector, parents=[latent_var], parametrization={"probs": CategoricalSelector(in_latent=self.latent_size, out_concepts=len(self.task_names), out_exogenous=self.n_rules, selector_hidden_layers=self.selector_hidden_layers)})
+
+        roles = EmbeddingVariable("rule_roles", distribution=OneHotCategorical, shape=(len(self.task_names), self.n_rules, n_concepts, 3))
+        roles_cpd = ParametricCPD(roles, parents=[], parametrization={"probs": RuleMemory(len(self.task_names), self.n_rules, n_concepts, self.memory_latent_size, self.memory_decoder_hidden_layers, hard_at_eval=self.hard_roles_at_eval)})
+
+        def aggregate(values):
+            return {"concepts": torch.cat([values[parent] for parent in list(concepts)], dim=-1), "selector": values[selector], "roles": values[roles]}
+
+        tasks = self.build_concept_variables(self.task_names, "tasks")
+        assert len(tasks) == 1, "CMR requires homogeneous binary task variables."
+        task_cpd = ParametricCPD(tasks[0], parents=[*concepts, selector, roles], parametrization={"probs": RuleTaskPredictor(out_concepts=tasks[0].size, in_concepts=n_concepts)}, aggregate=aggregate)
+
+        rec_tasks = EmbeddingVariable("tasks_with_rec", distribution=Bernoulli, shape=tasks[0].shape)
+        rec_cpd = ParametricCPD(rec_tasks, parents=[*concepts, selector, roles], parametrization={"probs": RuleReconstructionPredictor(out_concepts=tasks[0].size, in_concepts=n_concepts, rec_weight=self.rec_weight)}, aggregate=aggregate)
+
+        return BayesianNetwork(variables=[input_var, latent_var, *concepts, selector, roles, *tasks, rec_tasks], factors=[input_cpd, latent_cpd, *concept_cpds, selector_cpd, roles_cpd, task_cpd, rec_cpd])
diff --git a/torch_concepts/nn/modules/loss.py b/torch_concepts/nn/modules/loss.py
index 670d5593..fe13a412 100644
--- a/torch_concepts/nn/modules/loss.py
+++ b/torch_concepts/nn/modules/loss.py
@@ -597,4 +597,52 @@ def forward(
             mask = torch.isfinite(input)
         if mask.any():
             return self.scale * input[mask].abs().mean()
-        return torch.tensor(0.0, device=input.device)
\ No newline at end of file
+        return torch.tensor(0.0, device=input.device)
+
+
+class CMRBlendedLoss(TypeAwareLoss):
+    """CMR objective that switches task path based on the binary label.
+
+    Negative examples (``y=0``) supervise the ordinary task path ``y_pred``.
+    Positive examples (``y=1``) supervise the reconstruction-aware task path
+    ``y_pred_with_rec``.
+    Concept supervision remains standard BCE on the intermediate concepts.
+
+    The reconstruction-aware task prediction is expected in
+    ``output.extra["input_with_rec"]``.
+    """
+
+    def __init__(self, task_names, concept_weight: float = 1.0, task_weight: float = 1.0):
+        super().__init__()
+        self.task_names = list(task_names)
+        self.concept_weight = float(concept_weight)
+        self.task_weight = float(task_weight)
+
+    def forward(self, output: ModelOutput, target=None) -> torch.Tensor:
+        target = target if target is not None else output.target
+        if target is None:
+            raise ValueError("CMRBlendedLoss requires a concept-space target.")
+        if output.probs is None:
+            raise ValueError("CMRBlendedLoss requires Bernoulli probability outputs.")
+        if not output.extra or "task_input" not in output.extra or "input_with_rec" not in output.extra:
+            raise ValueError("CMRBlendedLoss requires output.extra['task_input'] and output.extra['input_with_rec'].")
+
+        task_target = target[self.task_names].to(output.probs.dtype)
+        task_pred = output.extra["task_input"].to(output.probs.dtype)
+        rec_pred = output.extra["input_with_rec"].to(task_pred.dtype)
+        if task_pred.shape != rec_pred.shape or task_pred.shape != task_target.shape:
+            raise ValueError("CMR task predictions and targets must have identical shapes.")
+
+        concept_names = [name for name in target.annotation.labels if name not in self.task_names]
+        if concept_names:
+            concept_loss = nn.functional.binary_cross_entropy(
+                output.probs[concept_names], target[concept_names].to(output.probs.dtype)
+            )
+        else:
+            concept_loss = task_pred.new_zeros(())
+
+        normal_bce = nn.functional.binary_cross_entropy(task_pred, task_target, reduction="none")
+        rec_bce = nn.functional.binary_cross_entropy(rec_pred, task_target, reduction="none")
+        switched = (1.0 - task_target) * normal_bce + task_target * rec_bce
+        task_loss = switched.mean()
+        return self.concept_weight * concept_loss + self.task_weight * task_loss
diff --git a/torch_concepts/nn/modules/low/encoders/selector.py b/torch_concepts/nn/modules/low/encoders/selector.py
new file mode 100644
index 00000000..bd374862
--- /dev/null
+++ b/torch_concepts/nn/modules/low/encoders/selector.py
@@ -0,0 +1,108 @@
+"""
+Memory selector module for memory selection.
+
+This module provides a memory-based selector that learns to attend over
+a memory bank of concept exogenous.
+"""
+import numpy as np
+import torch
+import torch.nn.functional as F
+
+
+from ..base.layer import BaseConceptLayer
+
+
+class CategoricalSelector(BaseConceptLayer):
+    """
+    Categorical selector that outputs concept-wise assignment probabilities.
+
+    This module maps latent inputs to logits of shape
+    ``(batch_size, out_concepts, out_exogenous)`` and applies a softmax over
+    the exogenous dimension to produce normalized mixing probabilities.
+
+    Attributes:
+        out_exogenous (int): Hidden width used in the selector MLP.
+        out_concepts (int): Number of output concepts.
+        selector_hidden_layers (int): Number of hidden layers in the selector MLP.
+        selector (nn.Sequential): Attention network for memory selection.
+
+    Args:
+        in_latent: Number of input latent features.
+        out_exogenous: Number of output exogenous features.
+        out_concepts: Number of output concept representations.
+        selector_hidden_layers: Number of hidden layers in the selector MLP.
+            Must be >= 0.
+        *args: Additional positional arguments for linear layers in the selector.
+        **kwargs: Additional keyword arguments for linear layers in the selector.
+
+    References:
+        Debot et al. "Interpretable Concept-Based Memory Reasoning", NeurIPS 2024. https://arxiv.org/abs/2407.15527
+    """
+    def __init__(
+        self,
+        in_latent: int,
+        out_exogenous: int,  # nb_rules
+        out_concepts: int,  # nb_tasks
+        selector_hidden_layers: int = 1,
+        *args,
+        **kwargs,
+    ):
+        """
+        Initialize the categorical selector.
+
+        Args:
+            in_latent: Number of input latent features.
+            out_exogenous: Number of output exogenous features.
+            out_concepts: Number of output concepts.
+            selector_hidden_layers: Number of hidden layers in the selector
+                MLP. Must be >= 0.
+            *args: Additional positional arguments for linear layers in the selector.
+            **kwargs: Additional keyword arguments for linear layers in the selector.
+        """
+        super().__init__(
+            in_embeddings=in_latent,
+            out_concepts=out_concepts,
+        )
+        self.in_latent = in_latent
+        if selector_hidden_layers < 0:
+            raise ValueError("selector_hidden_layers must be >= 0")
+
+        self.out_exogenous = out_exogenous
+        self.out_concepts = out_concepts
+        self.selector_hidden_layers = selector_hidden_layers
+        self._selector_out_shape = (out_concepts, out_exogenous)
+        self._selector_out_dim = np.prod(self._selector_out_shape).item()
+
+        selector_layers = []
+        in_features = in_latent
+        for _ in range(selector_hidden_layers):
+            selector_layers.extend([
+                torch.nn.Linear(in_features, in_latent, *args, **kwargs),
+                torch.nn.ReLU(),
+            ])
+            in_features = in_latent
+        selector_layers.extend([
+            torch.nn.Linear(in_features, self._selector_out_dim, *args, **kwargs),
+            torch.nn.Unflatten(-1, self._selector_out_shape),
+        ])
+        self.selector = torch.nn.Sequential(*selector_layers)
+
+    def forward(
+        self,
+        latent: torch.Tensor,
+    ) -> torch.Tensor:
+        """
+        Compute concept-wise mixing probabilities from latent input.
+
+        Applies the selector MLP and normalizes logits with softmax over
+        the exogenous axis.
+
+        Args:
+            latent: Input latent of shape (batch_size, in_latent).
+
+        Returns:
+            torch.Tensor: Mixing probabilities of shape
+                         (batch_size, out_concepts, out_exogenous).
+        """
+        mixing_coeff = self.selector(latent)
+        return torch.softmax(mixing_coeff, dim=-1)  # [Batch x Task x Memory]
diff --git a/torch_concepts/nn/modules/low/predictors/rule.py b/torch_concepts/nn/modules/low/predictors/rule.py
new file mode 100644
index 00000000..1e9ec06e
--- /dev/null
+++ b/torch_concepts/nn/modules/low/predictors/rule.py
@@ -0,0 +1,80 @@
+import torch
+
+from ..base.layer import BaseConceptLayer
+from ....functional import grouped_concept_exogenous_mixture, replace_expand_cols
+from typing import List
+
+
+class RuleMemory(torch.nn.Module):
+    """Learnable rule memory decoded into categorical role probabilities.
+
+    During training the decoded roles remain soft probabilities. During eval,
+    ``hard_at_eval=True`` converts each 3-way role categorical to its argmax
+    one-hot mode.
+
+    References:
+        Debot et al. "Interpretable Concept-Based Memory Reasoning", NeurIPS 2024.
+        https://arxiv.org/abs/2407.15527
+    """
+    def __init__(self, n_tasks, n_rules, n_concepts, latent_size=100, hidden_layers=1, hard_at_eval=False):
+        super().__init__()
+        self.hard_at_eval = hard_at_eval
+        self.shape = (n_tasks, n_rules, n_concepts, 3)
+        width = n_rules * n_concepts * 3
+        self.memory = torch.nn.Embedding(n_tasks, latent_size)
+        layers = [torch.nn.Linear(latent_size, width)]
+        for _ in range(hidden_layers):
+            layers += [torch.nn.LeakyReLU(), torch.nn.Linear(width, width)]
+        layers += [torch.nn.Unflatten(-1, (n_rules, n_concepts, 3))]
+        self.decoder = torch.nn.Sequential(*layers)
+    def forward(self):
+        pred = torch.softmax(self.decoder(self.memory.weight), dim=-1)
+        if (not self.training) and self.hard_at_eval:
+            idx = pred.argmax(dim=-1)
+            pred = torch.nn.functional.one_hot(idx, num_classes=pred.shape[-1]).to(pred.dtype)
+        assert torch.all((pred >= 0) & (pred <= 1)), "Decoded memory should be in [0, 1]"
+        return pred
+
+
+class RuleTaskPredictor(BaseConceptLayer):
+    """Compute the ordinary CMR task probability from concepts, selector and roles.
+
+    References:
+        Debot et al. "Interpretable Concept-Based Memory Reasoning", NeurIPS 2024.
+        https://arxiv.org/abs/2407.15527
+    """
+    def forward(self, concepts, selector, roles):
+        c = concepts.detach().unsqueeze(1).unsqueeze(1)
+        per_rule = (c * roles[..., 0] + (1.0 - c) * roles[..., 1] + roles[..., 2]).prod(dim=-1)
+        pred = (per_rule * selector).sum(dim=-1)
+        eps = 0.0001
+        pred = eps + (1 - 2 * eps) * pred  # numerical stability
+        return pred
+
+
+class RuleReconstructionPredictor(BaseConceptLayer):
+    """Compute the reconstruction-aware CMR task probability.
+
+    For each rule, this predictor multiplies its task satisfaction probability
+    by its reconstruction probability raised to ``rec_weight``. A weight of
+    zero disables reconstruction within this branch; larger non-negative
+    weights make reconstruction agreement more influential.
+
+    References:
+        Debot et al. "Interpretable Concept-Based Memory Reasoning", NeurIPS 2024.
+        https://arxiv.org/abs/2407.15527
+    """
+    def __init__(self, rec_weight=1.0, **kwargs):
+        super().__init__(**kwargs)
+        if rec_weight < 0:
+            raise ValueError("rec_weight must be non-negative.")
+        self.rec_weight = rec_weight
+    def forward(self, concepts, selector, roles):
+        c = concepts.detach().unsqueeze(1).unsqueeze(1)
+        task_per_rule = (c * roles[..., 0] + (1.0 - c) * roles[..., 1] + roles[..., 2]).prod(dim=-1)
+        reconstruction_per_rule = (c * roles[..., 0] + (1.0 - c) * roles[..., 1] + 0.5 * roles[..., 2]).prod(dim=-1)
+        reconstruction_per_rule = torch.pow(reconstruction_per_rule + 1e-6, self.rec_weight)
+        pred = (task_per_rule * reconstruction_per_rule * selector).sum(dim=-1)
+        eps = 0.0001
+        pred = eps + (1 - 2 * eps) * pred  # numerical stability
+        return pred
diff --git a/torch_concepts/nn/modules/mid/inference/__init__.py b/torch_concepts/nn/modules/mid/inference/__init__.py
index a17d6eeb..49e70f1d 100644
--- a/torch_concepts/nn/modules/mid/inference/__init__.py
+++ b/torch_concepts/nn/modules/mid/inference/__init__.py
@@ -3,4 +3,4 @@
 Backends live in the :mod:`.torch` (pure-PyTorch) and :mod:`.pyro` subpackages.
 The public API is re-exported from :mod:`torch_concepts.nn`.
 """
-__all__: list[str] = []
+__all__: list[str] = []
\ No newline at end of file