new startHEADmaster

Now I tested and discussed with claude, and had a working version of
this project.  The rest is to port the claude version to my own
writings.  So this serves as the starting point before that attempt.

In the new attempt the design is better, and a number of bugs were
found and supposedly fixed, so let us see how this goes at the end.

8 files changed, 2531 insertions, 393 deletions

diff --git a/forest/src/archive_default.txt b/forest/src/archive_default.txt
index 5572ecc..d6e4e2e 100644
--- a/forest/src/archive_default.txt
+++ b/forest/src/archive_default.txt
@@ -523,3 +523,84 @@ mod forest_test {
         Ok(())
     }
 }
+
+// #[derive(Debug, Clone, Copy)]
+// enum Item<T: Copy + Clone + std::fmt::Debug> {
+//     Seen(T),
+//     Unseen(T),
+// }
+
+// use Item::*;
+
+// impl<T: Copy + Clone + std::fmt::Debug> Item<T> {
+//     fn is_seen(self) -> bool {
+//         matches!(self, Self::Seen(_))
+//     }
+
+//     fn value(&self) -> T {
+//         match self {
+//             Seen(v) => *v,
+//             Unseen(v) => *v,
+//         }
+//     }
+// }
+
+// let query_len = self.query.nodes_len();
+
+// let forest_len = forest.nodes_len();
+
+// let mut stack: Vec<(SkipState, usize)> = Vec::with_capacity(forest_len);
+
+// stack.push((SkipState::Any, root));
+
+// let mut query_stack: Vec<usize> = Vec::with_capacity(query_len);
+
+// query_stack.push(query_root);
+
+// let mut match_result = false;
+
+// 'search: while let Some(item) = stack.pop() {
+//     // if item.is_seen() {
+//     //     if match_result {
+//     //         continue;
+//     //     }
+
+//     //     // fallback
+//     // }
+
+//     let (state, top) = item;
+
+//     let qtop = if let Some(qtop) = query_stack.last() {
+//         *qtop
+//     } else {
+//         panic!("Query stack becomes empty");
+//     };
+
+//     let top_label = forest.vertex_label(top)?.ok_or(QRError::NodeNoLabel(top))?;
+
+//     // let top_status = top_label.status();
+
+//     if top_label.is_packed() {
+//         self.packs.push((top, 0));
+
+//         stack.try_reserve(forest.degree(top)?)?;
+
+//         stack.extend(forest.children_of(top)?.map(|child| (state, child)));
+
+//         continue 'search;
+//     }
+
+//     let query_label = self.query.label_of_vertex(qtop)?;
+
+//     if query_label.matches(top_label.label(), top_label.status()) {
+//         // Success branch
+
+//         match_result = true;
+
+//         // TODO: Deal with captures later
+//     } else {
+//         // Failure branch
+
+//         match_result = false;
+//     }
+// }
diff --git a/forest/src/bytes.rs b/forest/src/bytes.rs
index d0d7e1e..6e854c8 100644
--- a/forest/src/bytes.rs
+++ b/forest/src/bytes.rs
@@ -65,18 +65,18 @@
 //! Each label has 3 parts:
 //!
 //! 1. Status: either Packed, Cloned, or Plain.  If the node is
-//! cloned, it has a clone index.  So in total this part occupies 1
-//! byte for the status and 8 bytes for the clone index.
+//!    cloned, it has a clone index.  So in total this part occupies 1
+//!    byte for the status and 8 bytes for the clone index.
 //! 2. Start and end: the range in the input sentence.  We just store
-//! two numbers here.  Hence this part occupies 16 bytes.
+//!    two numbers here.  Hence this part occupies 16 bytes.
 //! 3. Grammar label: either a terminal, a non-terminal, or a rule.
-//! Each variant needs a number to speify its index.  So in total this
-//! part occupies 1 byte for the variant discriminant and 8 bytes for
-//! the number.
+//!    Each variant needs a number to speify its index.  So in total
+//!    this part occupies 1 byte for the variant discriminant and 8
+//!    bytes for the number.
 //!
 //! To sum up, each label occupies 34 bytes.
 
-use super::Forest;
+use super::{Forest, ForestLabel, ForestLabelType};
 use grammar::{GrammarLabel, GrammarLabelType, TNT};
 use graph::error::Error as GError;
 
@@ -95,13 +95,16 @@ pub fn forest_to_bytes<F: Forest<GrammarLabel>>(forest: &F) -> Vec<u8> {
 
     let total_degree: usize = degrees.iter().copied().sum();
 
-    let len: usize = 8          // total number of bytes at the start
+    let label_offset: usize = 8 // total number of bytes = len
         + 3                     // special mark
         + 8                     // number of nodes
         + 8                     // offset of labels
         + 16 * nodes_len        // degree & offset for each node
         + 8 * total_degree      // children of each node
-        + 34 * nodes_len        // labels
+        ;
+
+    let len: usize = label_offset // see that variable for details.
+        + 34 * nodes_len          // labels
         ;
 
     // Then fill in the bytes.
@@ -113,17 +116,18 @@ pub fn forest_to_bytes<F: Forest<GrammarLabel>>(forest: &F) -> Vec<u8> {
     bytes.extend(len.to_be_bytes());
     bytes.extend([114, 101, 112]); // rep
     bytes.extend(nodes_len.to_be_bytes());
-    bytes.extend((len - 34 * nodes_len).to_be_bytes());
+    bytes.extend(label_offset.to_be_bytes());
 
-    let mut accumulated: usize = 0;
+    let mut accumulated: usize = 8 + 3 + 8 + 8 + 16 * nodes_len;
 
     // Then degrees and offsets
 
     for degree in degrees.iter().copied() {
         bytes.extend(degree.to_be_bytes());
-        bytes.extend((8 + 3 + 8 + 8 + 16 * nodes_len + 8 * accumulated).to_be_bytes());
+        bytes.extend(accumulated.to_be_bytes());
+        // bytes.extend((8 + 3 + 8 + 8 + 16 * nodes_len + 8 * accumulated).to_be_bytes());
 
-        accumulated += degree;
+        accumulated += 8 * degree;
     }
 
     // Then the children
@@ -337,43 +341,147 @@ impl TryFrom<Vec<u8>> for BytesForest {
 pub struct BytesChildrenIter<'a> {
     vec: &'a [u8],
     offset: usize,
+    label_offset: usize,
     degree: usize,
     current: usize,
 }
 
+/// The delegation of iterators for children and labels of nodes.
+#[derive(Debug, Clone)]
+pub struct BytesLabelChildrenIter<'a>(BytesChildrenIter<'a>);
+
+impl<'a> Iterator for BytesLabelChildrenIter<'a> {
+    type Item = (ForestLabel<GrammarLabel>, usize);
+
+    #[inline]
+    fn next(&mut self) -> Option<Self::Item> {
+        let iter = &self.0;
+
+        if iter.current >= iter.degree {
+            return None;
+        }
+
+        let len = iter.vec.len();
+
+        if iter.offset + 8 * iter.degree >= len {
+            dbg!();
+            return None;
+        }
+
+        if iter.label_offset + 34 * iter.degree > len {
+            dbg!();
+            return None;
+        }
+
+        let result_index = usize::from_be_bytes(
+            (&iter.vec[(iter.offset + 8 * iter.current)..(iter.offset + 8 * iter.current + 8)])
+                .try_into()
+                .ok()?,
+        );
+
+        let result_label = read_label(
+            &iter.vec[(iter.label_offset + iter.current * 34)
+                ..(iter.label_offset + iter.current * 34 + 34)],
+        )
+        .ok()?;
+
+        self.0.current += 1;
+
+        Some((result_label, result_index))
+    }
+
+    #[inline]
+    fn size_hint(&self) -> (usize, Option<usize>) {
+        (self.0.degree, Some(self.0.degree))
+    }
+}
+
+impl<'a> ExactSizeIterator for BytesLabelChildrenIter<'a> {
+    #[inline]
+    fn len(&self) -> usize {
+        self.0.degree
+    }
+}
+
 impl<'a> Iterator for BytesChildrenIter<'a> {
     type Item = usize;
 
+    #[inline]
     fn next(&mut self) -> Option<Self::Item> {
-        if self.current + 1 >= self.degree {
+        if self.current >= self.degree {
             return None;
         }
 
         let len = self.vec.len();
 
-        if self.offset + 8 * self.degree > len {
+        if self.offset + 8 * self.degree >= len {
             dbg!();
             return None;
         }
 
-        Some(usize::from_be_bytes(
+        let result = usize::from_be_bytes(
             (&self.vec[(self.offset + 8 * self.current)..(self.offset + 8 * self.current + 8)])
                 .try_into()
                 .ok()?,
-        ))
+        );
+
+        self.current += 1;
+
+        Some(result)
     }
 
+    #[inline]
     fn size_hint(&self) -> (usize, Option<usize>) {
         (self.degree, Some(self.degree))
     }
 }
 
 impl<'a> ExactSizeIterator for BytesChildrenIter<'a> {
+    #[inline]
     fn len(&self) -> usize {
         self.degree
     }
 }
 
+fn read_label(data: &[u8]) -> Result<ForestLabel<GrammarLabel>, GError> {
+    let len = data.len();
+
+    if len < 34 {
+        return Err(GError::IndexOutOfBounds(34, len));
+    }
+
+    let pack_status: ForestLabelType = match data.first().copied().unwrap() {
+        0 => ForestLabelType::Plain,
+        1 => ForestLabelType::Packed,
+        _ => ForestLabelType::Cloned(usize::from_be_bytes(
+            data[1..9].try_into().expect("should not happen"),
+        )),
+    };
+
+    let start = usize::from_be_bytes(data[9..17].try_into().expect("should not happen"));
+    let end = usize::from_be_bytes(data[17..25].try_into().expect("should not happen"));
+
+    let grammarlabeltype = match data.get(25).unwrap() {
+        0 => GrammarLabelType::TNT(TNT::Ter(usize::from_be_bytes(
+            data[26..34].try_into().expect("should not happen"),
+        ))),
+        1 => GrammarLabelType::TNT(TNT::Non(usize::from_be_bytes(
+            data[26..34].try_into().expect("should not happen"),
+        ))),
+        _ => GrammarLabelType::Rule(usize::from_be_bytes(
+            data[26..34].try_into().expect("should not happen"),
+        )),
+    };
+
+    let label = if end == 0 {
+        GrammarLabel::new(grammarlabeltype, start)
+    } else {
+        GrammarLabel::new_closed(grammarlabeltype, start, end)
+    };
+
+    Ok(ForestLabel::new(label, pack_status))
+}
+
 impl graph::Graph for BytesForest {
     type Iter<'a> = BytesChildrenIter<'a>
     where
@@ -404,6 +512,8 @@ impl graph::Graph for BytesForest {
                 .unwrap(),
         );
 
+        let label_offset = usize::from_be_bytes((&self.bytes[19..27]).try_into().unwrap());
+
         let current = 0;
 
         let vec = &self.bytes;
@@ -411,6 +521,7 @@ impl graph::Graph for BytesForest {
         Ok(BytesChildrenIter {
             vec,
             offset,
+            label_offset,
             degree,
             current,
         })
@@ -453,4 +564,118 @@ impl graph::Graph for BytesForest {
     }
 }
 
-// TODO: Implement labelgraph
+impl graph::VLabelGraph<ForestLabel<GrammarLabel>> for BytesForest {
+    type TupleIter<'a> = BytesLabelChildrenIter<'a>
+    where
+        Self: 'a,
+        ForestLabel<GrammarLabel>: 'a;
+
+    fn verlabren_of(&self, node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
+        if node_id >= self.nodes_len {
+            return Err(GError::IndexOutOfBounds(node_id, self.nodes_len));
+        }
+
+        let degree = usize::from_be_bytes(
+            (&self.bytes[(27 + 16 * node_id)..(35 + 16 * node_id)])
+                .try_into()
+                .unwrap(),
+        );
+
+        let offset = usize::from_be_bytes(
+            (&self.bytes[(35 + 16 * node_id)..(43 + 16 * node_id)])
+                .try_into()
+                .unwrap(),
+        );
+
+        let label_offset = usize::from_be_bytes((&self.bytes[19..27]).try_into().unwrap());
+
+        let current = 0;
+
+        let vec = &self.bytes;
+
+        Ok(BytesLabelChildrenIter(BytesChildrenIter {
+            vec,
+            offset,
+            label_offset,
+            degree,
+            current,
+        }))
+    }
+}
+
+impl graph::LabelGraph<ForestLabel<GrammarLabel>> for BytesForest {
+    fn query_label(&self, label: ForestLabel<GrammarLabel>) -> Option<usize> {
+        use graph::Graph;
+        let nodes_len = self.nodes_len();
+
+        for index in 0..nodes_len {
+            let index_label = read_label(
+                &self.bytes[(self.offset + index * 34)..(self.offset + index * 34 + 34)],
+            )
+            .ok()?;
+
+            if index_label == label {
+                return Some(index);
+            }
+        }
+
+        None
+    }
+
+    fn vertex_label(&self, node_id: usize) -> Result<Option<ForestLabel<GrammarLabel>>, GError> {
+        read_label(&self.bytes[(self.offset + node_id * 34)..(self.offset + node_id * 34 + 34)])
+            .map(Some)
+    }
+
+    type Iter<'a> = std::iter::Empty<usize>
+    where
+        Self: 'a;
+
+    type TupleIter<'a>=std::iter::Empty<(ForestLabel<GrammarLabel>, usize)>
+    where
+        Self: 'a,
+        ForestLabel<GrammarLabel>: 'a;
+
+    type LabelIter<'a> = std::iter::Empty<(&'a ForestLabel<GrammarLabel>, std::iter::Empty<usize>)>
+    where
+        Self: 'a,
+        ForestLabel<GrammarLabel>: 'a;
+
+    type EdgeLabelIter<'a> = std::iter::Empty<ForestLabel<GrammarLabel>>
+    where
+        Self: 'a,
+        ForestLabel<GrammarLabel>: 'a;
+
+    fn edge_label(
+        &self,
+        _source: usize,
+        _target: usize,
+    ) -> Result<Self::EdgeLabelIter<'_>, GError> {
+        unimplemented!("bytes forests do not have edge labels")
+    }
+
+    fn find_children_with_label(
+        &self,
+        _node_id: usize,
+        _label: &ForestLabel<GrammarLabel>,
+    ) -> Result<<Self as graph::LabelGraph<ForestLabel<GrammarLabel>>>::Iter<'_>, GError> {
+        unimplemented!("bytes forests do not have edge labels")
+    }
+
+    fn labels_of(&self, _node_id: usize) -> Result<Self::LabelIter<'_>, GError> {
+        unimplemented!("bytes forests do not have edge labels")
+    }
+
+    fn tuples_of(&self, _node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
+        unimplemented!("bytes forests do not have edge labels")
+    }
+
+    fn has_edge_label(
+        &self,
+        _node_id: usize,
+        _label: &ForestLabel<GrammarLabel>,
+        _target: usize,
+    ) -> Result<bool, GError> {
+        unimplemented!("bytes forests do not have edge labels")
+    }
+}
diff --git a/forest/src/iter/archive.txt b/forest/src/iter/archive.txt
new file mode 100644
index 0000000..dcdae7f
--- /dev/null
+++ b/forest/src/iter/archive.txt
@@ -0,0 +1,991 @@
+/// Archive of failed attempts
+
+/// The label of a query graph.
+///
+/// Such a label is responsible for matching forest labels.
+pub trait QueryLabel: GraphLabel {
+    /// The type of labels this matches against.
+    type Obj: GraphLabel;
+
+    /// Determine whether or not a label matches this query label.
+    fn matches(&self, label: Self::Obj, status: ForestLabelType) -> bool {
+        let _ = label;
+        let _ = status;
+
+        false
+    }
+
+    /// Perform the query for one single query label.
+    fn run<F>(
+        &self,
+        forest: &F,
+        node: usize,
+        packs: &mut Vec<(usize, usize)>,
+    ) -> Result<Option<usize>, QRError>
+    where
+        F: LabelGraph<ForestLabel<<Self as QueryLabel>::Obj>>,
+    {
+        let mut dfs_stack: Vec<usize> = vec![node];
+
+        let nodes_len = forest.nodes_len();
+
+        while let Some(node) = dfs_stack.pop() {
+            let label = match forest.vertex_label(node) {
+                Ok(None) => return Err(QRError::NodeNoLabel(node)),
+                Err(_) => return Err(QRError::IndexOutOfBounds(node, nodes_len)),
+                Ok(Some(label)) => label,
+            };
+
+            let status = label.status();
+            let label = label.label();
+
+            if status == ForestLabelType::Packed {
+                packs.push((node, 0));
+
+                dfs_stack.extend(forest.children_of(node).unwrap().take(1));
+
+                continue;
+            }
+
+            if self.matches(label, status) {
+                return Ok(Some(node));
+            }
+
+            dfs_stack.extend(
+                forest
+                    .children_of(node)
+                    .unwrap()
+                    .collect::<Vec<_>>()
+                    .into_iter()
+                    .rev(),
+            );
+        }
+
+        Ok(None)
+    }
+}
+
+/// How should the current query skip nodes
+#[derive(Debug, Clone, Copy, Default, Ord, PartialOrd, Hash, Eq, PartialEq)]
+pub enum SkipState {
+    /// Skip anything in the depth-first search
+    Any,
+    /// Skip siblings
+    Sibling,
+    /// Do not skip
+    #[default]
+    None,
+}
+
+impl std::fmt::Display for SkipState {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::Any => write!(f, "skip anything"),
+            Self::Sibling => write!(f, "skip siblings"),
+            Self::None => write!(f, "no skipping"),
+        }
+    }
+}
+
+impl From<bool> for SkipState {
+    fn from(value: bool) -> Self {
+        if value {
+            Self::Sibling
+        } else {
+            Self::None
+        }
+    }
+}
+
+/// A branch switch is of the form (n, m, v), where n means the n-th
+/// node in the forest and m means the m-th child, and v is the vector
+/// of children of the packed node.
+#[derive(Debug, Clone, Default)]
+pub struct BranchSwitch {
+    n: usize,
+    m: usize,
+    v: Vec<usize>,
+}
+
+impl BranchSwitch {
+    #[allow(dead_code)]
+    fn new(n: usize, m: usize, v: Vec<usize>) -> Self {
+        Self { n, m, v }
+    }
+
+    #[allow(dead_code)]
+    fn node(&self) -> usize {
+        self.n
+    }
+
+    #[allow(dead_code)]
+    fn child(&self) -> Option<usize> {
+        self.v.get(self.m).copied()
+    }
+
+    #[allow(dead_code)]
+    fn advance(&mut self) -> bool {
+        if self.m < self.v.len() {
+            self.m += 1;
+
+            return true;
+        }
+
+        false
+    }
+}
+
+// REVIEW: Change the types into boxed slices so that we can guarantee
+// no capacity is ever wasted.
+
+/// The type of forests that would be walked.
+///
+/// In plain terms this is a preprocessed graph that provides some
+/// extra information than a plain graph, which will be useful in
+/// walking.
+#[derive(Debug, Clone)]
+#[allow(unused)]
+pub struct QForest {
+    current: Option<usize>,
+    /// Just for example uses.
+    pub dfs_next: Vec<usize>,
+    /// Just for example uses.
+    pub subtree_ranges: Vec<std::ops::RangeInclusive<usize>>,
+    // /// Just for example uses.
+    // pub next_branch: Vec<Option<usize>>,
+    num_seen_switches: usize,
+    /// Switches of branches are packed nodes along with the current
+    /// index of the child.  One can switch to use a different child
+    /// to "switch branches".
+    pub branch_switches: Vec<BranchSwitch>,
+    /// Map ordinary node numbering to depth-first numbering.
+    ord_to_dfs: Vec<usize>,
+    /// Map depth-first numbering to ordinary node numbering.
+    dfs_to_ord: Vec<usize>,
+    restriction: Option<usize>,
+}
+
+impl QForest {
+    fn new(
+        dfs_next: Vec<usize>,
+        subtree_ranges: Vec<std::ops::RangeInclusive<usize>>,
+        // next_branch: Vec<Option<usize>>,
+        branch_switches: Vec<BranchSwitch>,
+        ord_to_dfs: Vec<usize>,
+        dfs_to_ord: Vec<usize>,
+    ) -> Self {
+        let current: Option<usize> = None;
+        let restriction: Option<usize> = None;
+        let num_seen_switches = 0usize;
+
+        Self {
+            current,
+            dfs_next,
+            subtree_ranges,
+            // next_branch,
+            num_seen_switches,
+            branch_switches,
+            restriction,
+            ord_to_dfs,
+            dfs_to_ord,
+        }
+    }
+
+    /// Translate the node index from ordinary index to depth-first
+    /// index.
+    pub fn translate_to_dfs(&self, id: usize) -> Option<usize> {
+        self.ord_to_dfs.get(id).copied()
+    }
+
+    /// Translate the node index from depth-first index to ordinary
+    /// index.
+    pub fn translate_from_dfs(&self, id: usize) -> Option<usize> {
+        self.dfs_to_ord.get(id).copied()
+    }
+}
+
+impl Default for QForest {
+    fn default() -> Self {
+        let current = Default::default();
+        let dfs_next = Default::default();
+        let subtree_ranges = Default::default();
+        // let next_branch = Default::default();
+        let branch_switches = Default::default();
+        let restriction = Default::default();
+        let num_seen_switches = 0usize;
+
+        let ord_to_dfs = Default::default();
+        let dfs_to_ord = Default::default();
+
+        Self {
+            current,
+            dfs_next,
+            subtree_ranges,
+            // next_branch,
+            num_seen_switches,
+            branch_switches,
+            restriction,
+            ord_to_dfs,
+            dfs_to_ord,
+        }
+    }
+}
+
+impl Queryable for QForest {
+    fn current(&self) -> Option<usize> {
+        self.current
+    }
+
+    fn next(&mut self) -> bool {
+        let current = if let Some(current) = self.current {
+            current
+        } else {
+            return false;
+        };
+
+        let mut next = if let Some(next) = self.dfs_next.get(current) {
+            *next
+        } else {
+            return false;
+        };
+
+        match self.branch_switches.get(self.num_seen_switches) {
+            Some(switch) if switch.node() == next => {
+                if let Some(child) = switch.child() {
+                    next = child;
+                }
+
+                self.num_seen_switches += 1;
+            }
+            _ => {}
+        }
+
+        let dfs_len = self.dfs_next.len();
+
+        if (0..dfs_len).contains(&next) {
+            self.current.replace(next);
+
+            true
+        } else {
+            false
+        }
+    }
+
+    fn next_branch(&mut self) -> bool {
+        self.branch_switches.iter_mut().any(BranchSwitch::advance)
+
+        // let current = if let Some(current) = self.current {
+        //     current
+        // } else {
+        //     return false;
+        // };
+
+        // let next = if let Some(next) = self.next_branch.get(current).copied().flatten() {
+        //     next
+        // } else {
+        //     return false;
+        // };
+
+        // self.current.replace(next);
+
+        // true
+    }
+
+    fn restrict(&mut self, node: usize) -> bool {
+        let range = if let Some(range) = self.subtree_ranges.get(node) {
+            range
+        } else {
+            return false;
+        };
+
+        self.restriction.replace(node);
+
+        let current = if let Some(current) = &self.current {
+            current
+        } else {
+            return false;
+        };
+
+        range.contains(current)
+    }
+}
+
+/// Calculate some extra information for later processing.
+pub fn preprocess_forest<T, F>(forest: F, root: usize) -> Result<QForest, GError>
+where
+    T: GraphLabel,
+    F: LabelGraph<ForestLabel<T>>,
+{
+    let forest_len = forest.nodes_len();
+
+    // let current: Option<usize> = None;
+
+    let mut dfs_next: Vec<usize> = (0..forest_len).map(|i| i + 1).collect();
+
+    let mut subtree_ranges: Vec<std::ops::RangeInclusive<usize>> =
+        (0..forest_len).map(|i| i..=i).collect();
+
+    // let next_branch: Vec<Option<usize>> = std::iter::repeat(None).take(forest_len).collect();
+
+    let mut branch_switches: Vec<BranchSwitch> = Vec::with_capacity(forest_len);
+
+    // let restriction: Option<usize> = None;
+
+    let mut dfs_to_ord: Vec<usize> = (0..forest_len).collect();
+    let mut ord_to_dfs: Vec<usize> = (0..forest_len).collect();
+
+    if forest.is_empty_node(root).unwrap_or(true) {
+        if let Some(root_item) = dfs_next.get_mut(root) {
+            *root_item = forest_len; // out of bounds
+        }
+
+        return Ok(QForest::new(
+            // current,
+            dfs_next,
+            subtree_ranges,
+            // next_branch,
+            // restriction,
+            branch_switches,
+            ord_to_dfs,
+            dfs_to_ord,
+        ));
+    }
+
+    #[derive(Debug, Clone, Copy)]
+    enum StackItem {
+        Normal(usize),
+        Reflection(usize),
+    }
+
+    impl From<usize> for StackItem {
+        fn from(value: usize) -> Self {
+            Self::Normal(value)
+        }
+    }
+
+    impl StackItem {
+        fn is_normal(&self) -> bool {
+            matches!(self, Self::Normal(_))
+        }
+
+        fn get_value(&self) -> usize {
+            match *self {
+                StackItem::Normal(item) => item,
+                StackItem::Reflection(item) => item,
+            }
+        }
+
+        fn reflect(&self) -> Self {
+            match self {
+                StackItem::Normal(item) => Self::Reflection(*item),
+                StackItem::Reflection(_) => *self,
+            }
+        }
+    }
+
+    let mut stack: Vec<StackItem> = Vec::with_capacity(forest_len);
+
+    let mut prev = root;
+
+    let mut children_vec: Vec<usize> = Vec::with_capacity(forest_len);
+
+    stack.extend({
+        children_vec.extend(forest.children_of(root)?);
+
+        children_vec.iter().copied().rev().map(StackItem::from)
+    });
+
+    let mut seen_nodes: HashMap<usize, Option<usize>> = HashMap::with_capacity(forest_len);
+
+    seen_nodes.insert(root, None);
+
+    let mut dfs_order = 1usize;
+
+    while let Some(top) = stack.pop() {
+        let top_value = top.get_value();
+        let top_label = forest
+            .vertex_label(top_value)?
+            .expect("node without label: {top_value}");
+        if top.is_normal() {
+            if let Some(prev_item) = dfs_next.get_mut(prev) {
+                *prev_item = top_value;
+            }
+
+            if let std::collections::hash_map::Entry::Vacant(e) = seen_nodes.entry(top_value) {
+                if let Some(item) = dfs_to_ord.get_mut(dfs_order) {
+                    *item = top_value;
+                }
+
+                if let Some(item) = ord_to_dfs.get_mut(top_value) {
+                    *item = dfs_order;
+                }
+
+                dfs_order += 1;
+
+                stack.push(top.reflect());
+
+                e.insert(top_label.clone_index().map(|_| prev));
+
+                children_vec.clear();
+
+                children_vec.extend(forest.children_of(top_value)?);
+
+                stack.extend(children_vec.iter().copied().rev().map(StackItem::from));
+            }
+
+            prev = top_value;
+        } else {
+            if let Some(range) = subtree_ranges.get_mut(top_value) {
+                let top_dfs = ord_to_dfs.get(top_value).copied().unwrap();
+                let prev_dfs = ord_to_dfs.get(prev).copied().unwrap();
+
+                *range = top_dfs..=prev_dfs;
+            }
+
+            // if top_is_cloned {
+            //     if let Some(next) = stack.iter().rev().copied().find(StackItem::is_normal) {
+            //         // update branches
+            //         let top_dfs = ord_to_dfs.get(top_value).copied().unwrap();
+            //         let prev_dfs = ord_to_dfs.get(prev).copied().unwrap();
+
+            //         for i in (top_dfs + 1)..=prev_dfs {
+            //             let j = dfs_to_ord.get(i).copied().unwrap();
+
+            //             match next_branch.get_mut(j) {
+            //                 Some(branch) if branch.is_none() => {
+            //                     *branch = Some(next.get_value());
+            //                 }
+
+            //                 _ => {}
+            //             }
+            //         }
+            //     }
+            // }
+
+            if !top_label.is_packed() {
+                continue;
+            }
+
+            let top_children: Vec<usize> = forest.children_of(top_value)?.collect();
+
+            if let Some(next) = stack.iter().rev().copied().find(StackItem::is_normal) {
+                for child in top_children.iter() {
+                    let child_prev = if let Some(Some(item)) = seen_nodes.get(child) {
+                        *item
+                    } else {
+                        continue;
+                    };
+
+                    if let Some(prev_item) = dfs_next.get_mut(child_prev) {
+                        *prev_item = next.get_value();
+
+                        // dbg!(child_prev, next.get_value(), child);
+                    }
+                }
+            }
+
+            branch_switches.push(BranchSwitch::new(top_value, 0, top_children));
+        }
+    }
+
+    Ok(QForest::new(
+        // current,
+        dfs_next,
+        subtree_ranges,
+        // next_branch,
+        // restriction,
+        branch_switches,
+        ord_to_dfs,
+        dfs_to_ord,
+    ))
+}
+
+// TODO: The query graph needs the ability to jump to a parent and
+// restore the querying contexts.
+
+/// The type that represents queries as dependency graphs of query
+/// labels.
+///
+/// The labels of vertices are query labels that match against forest
+/// labels, whereas the labels of graphs indicate whether or not the
+/// edge should allow skipping siblings before matching.
+#[derive(Debug, Clone, Graph)]
+pub struct Query<T: QueryLabel> {
+    graph: SLGraph<SkipState>,
+    vertex_labels: Vec<T>,
+    root: Option<usize>,
+}
+
+impl<T: QueryLabel> Query<T> {
+    /// Construct a new query from raw fields.
+    pub fn new(graph: SLGraph<SkipState>, vertex_labels: Vec<T>, root: Option<usize>) -> Self {
+        Self {
+            graph,
+            vertex_labels,
+            root,
+        }
+    }
+
+    /// Return the root of this query.
+    ///
+    /// This returns `None` if and only if the graph is empty.
+    pub fn root(&self) -> Option<usize> {
+        self.root
+    }
+
+    /// Return the label of the given node.
+    ///
+    /// This cannot be implemented as a method of the `LabelGraph`
+    /// trait as the type differs from that of the labels.
+    #[inline]
+    pub fn label_of_vertex(&self, node: usize) -> Result<T, GError> {
+        self.vertex_labels
+            .get(node)
+            .copied()
+            .ok_or(GError::IndexOutOfBounds(node, self.vertex_labels.len()))
+    }
+}
+
+impl<T: QueryLabel> Default for Query<T> {
+    fn default() -> Self {
+        let graph = Default::default();
+        let vertex_labels = Vec::new();
+        let root = None;
+
+        Self {
+            graph,
+            vertex_labels,
+            root,
+        }
+    }
+}
+
+// delegation of labelgraph
+
+impl<T: QueryLabel> LabelGraph<SkipState> for Query<T> {
+    type Iter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::Iter<'a>
+    where
+        Self: 'a;
+
+    type TupleIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::TupleIter<'a>
+    where
+        Self: 'a;
+
+    type LabelIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::LabelIter<'a>
+    where
+        Self: 'a;
+
+    type EdgeLabelIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::EdgeLabelIter<'a>
+    where
+        Self: 'a;
+
+    #[inline]
+    fn edge_label(&self, source: usize, target: usize) -> Result<Self::EdgeLabelIter<'_>, GError> {
+        self.graph.edge_label(source, target)
+    }
+
+    fn find_children_with_label(
+        &self,
+        _node_id: usize,
+        _label: &SkipState,
+    ) -> Result<<Self as LabelGraph<SkipState>>::Iter<'_>, GError> {
+        unimplemented!("Unsupported operation")
+    }
+
+    fn labels_of(&self, _node_id: usize) -> Result<Self::LabelIter<'_>, GError> {
+        unimplemented!("Unsupported operation")
+    }
+
+    #[inline]
+    fn tuples_of(&self, node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
+        self.graph.tuples_of(node_id)
+    }
+
+    #[inline]
+    fn has_edge_label(
+        &self,
+        node_id: usize,
+        label: &SkipState,
+        target: usize,
+    ) -> Result<bool, GError> {
+        self.graph.has_edge_label(node_id, label, target)
+    }
+}
+
+#[derive(Debug, Clone)]
+/// The type of errors that can occur when querying forests.
+pub enum QRError {
+    /// An index is out of bounds.
+    IndexOutOfBounds(usize, usize),
+    /// A node has no labels while it is required to have one.
+    NodeNoLabel(usize),
+    /// Encounter an invalid error in converting from an error of
+    /// graphs.
+    InvalidGraphError(GError),
+    /// Not enough memory for vectors.
+    NoMemory(std::collections::TryReserveError),
+}
+
+impl std::fmt::Display for QRError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            QRError::IndexOutOfBounds(n, b) => {
+                write!(f, "index {n} out of bound {b}")
+            }
+            QRError::NodeNoLabel(n) => write!(f, "the node {n} has no label but should have one"),
+            QRError::InvalidGraphError(g) => {
+                write!(f, "an invalid error for forest iterators: {g}")
+            }
+            QRError::NoMemory(e) => write!(f, "not enough memory: {e}"),
+        }
+    }
+}
+
+impl std::error::Error for QRError {}
+
+impl From<GError> for QRError {
+    fn from(value: GError) -> Self {
+        match value {
+            GError::IndexOutOfBounds(n, b) => Self::IndexOutOfBounds(n, b),
+            _ => Self::InvalidGraphError(value),
+        }
+    }
+}
+
+impl From<std::collections::TryReserveError> for QRError {
+    fn from(value: std::collections::TryReserveError) -> Self {
+        Self::NoMemory(value)
+    }
+}
+
+/// An implementation of the forest iterator for querying purposes.
+#[derive(Debug, Clone)]
+pub struct QRun<T: QueryLabel> {
+    /// The currently visited node.
+    current: Option<Vec<usize>>,
+    /// The queries dependency graph.
+    query: Query<T>,
+    /// The nodes of the query graph, in the depth-first order.
+    query_dfs: Vec<usize>,
+    /// The nodes of the forest, in the depth-first order.
+    forest_dfs: Vec<usize>,
+    /// The packed nodes that can be used to change to another branch
+    /// if needed.
+    ///
+    /// The tuple is in the following format
+    ///
+    /// (forest_index, query_index),
+    ///
+    /// where the indicies are the indices of nodes, in the
+    /// depth-first order, of the respective graphs.
+    packs: Vec<(usize, usize)>,
+}
+
+impl<T: QueryLabel> Default for QRun<T> {
+    fn default() -> Self {
+        let current = Default::default();
+        let query = Default::default();
+        let packs = Default::default();
+
+        let query_dfs = Default::default();
+        let forest_dfs = Default::default();
+
+        Self {
+            current,
+            query,
+            query_dfs,
+            forest_dfs,
+            packs,
+        }
+    }
+}
+
+impl<T: QueryLabel> ForestIter for QRun<T> {
+    #[inline]
+    fn current(&self) -> Option<&[usize]> {
+        self.current.as_ref().map(|v| &v[..])
+    }
+
+    type ForestLabel = <T as QueryLabel>::Obj;
+
+    type Item<'a> = Option<&'a [usize]>
+    where
+        Self: 'a;
+
+    type Error = QRError;
+
+    #[allow(unused)]
+    #[allow(clippy::collapsible_if)]
+    fn next<F>(&mut self, forest: &F, root: usize) -> Result<State, Self::Error>
+    where
+        F: LabelGraph<ForestLabel<Self::ForestLabel>>,
+    {
+        #[derive(Debug, Clone, Copy, Default)]
+        enum QueryIndex {
+            /// A normal index
+            Normal(usize),
+            #[default]
+            /// This means the query is now accepting anything.
+            Empty,
+        }
+
+        macro_rules! query_match (
+            ($index:expr, $label:expr, $status:expr) => {
+                if let QueryIndex::Normal($index) = index {
+                    self.query
+                        .label_of_vertex(qtop)
+                        .map(|q| q.matches($label, $status))
+                } else {
+                    Ok(true)
+                }
+            }
+        );
+
+        #[derive(Debug, Clone, Copy)]
+        enum ForestIndex {
+            /// A normal index
+            Normal(usize),
+            /// How many times we fallback in the query graph
+            FallbackTimes(usize),
+        }
+
+        impl ForestIndex {
+            fn is_normal(&self) -> bool {
+                matches!(self, Self::Normal(_))
+            }
+
+            // /// Increment the index if it is normal, decrement it otherwise.
+            // ///
+            // /// Return `None` if the advanced index would be out of bounds.
+            // fn advance_forest_index(&self, forest_len: usize) -> Option<Self> {
+            //     match self {
+            //         ForestIndex::Normal(index) => {
+            //             if *index + 1 >= forest_len {
+            //                 None
+            //             } else {
+            //                 Some(Self::Normal(*index + 1))
+            //             }
+            //         }
+            //         ForestIndex::FallbackTimes(time) => {
+            //             if *time == 0 {
+            //                 None
+            //             } else {
+            //                 Some(Self::FallbackTimes(*time - 1))
+            //             }
+            //         }
+            //     }
+            // }
+
+            fn get_index(&self, current_indices: &[usize]) -> Result<usize, QRError> {
+                match self {
+                    ForestIndex::Normal(index) => Ok(*index),
+                    ForestIndex::FallbackTimes(time) => current_indices
+                        .get(*time)
+                        .copied()
+                        .ok_or(QRError::IndexOutOfBounds(*time, current_indices.len())),
+                }
+            }
+
+            fn get_query_index(&self) -> usize {
+                match self {
+                    ForestIndex::Normal(_) => 0,
+                    ForestIndex::FallbackTimes(time) => *time,
+                }
+            }
+        }
+
+        let query_empty_p = self.query.is_empty();
+
+        let query_len = self.query.nodes_len();
+        let forest_len = forest.nodes_len();
+
+        let query_root = if let Some(r) = self.query.root() {
+            r
+        } else {
+            return Ok(State::Finish);
+        };
+
+        let orig_top: ForestIndex = match &self.current {
+            Some(current) => {
+                // re-starting search
+                ForestIndex::FallbackTimes(current.len() - 1)
+            }
+            None => {
+                // fresh new search
+                ForestIndex::Normal(0)
+            }
+        };
+
+        let mut top = orig_top.get_index(if let Some(current) = &self.current {
+            &current[..]
+        } else {
+            &[]
+        })?;
+
+        let mut qtop = orig_top.get_query_index();
+
+        let mut candidate_captures: Vec<usize> = self
+            .current
+            .clone()
+            .unwrap_or_else(|| std::iter::repeat(0).take(query_len).collect());
+
+        'search: loop {
+            let top_index = self
+                .forest_dfs
+                .get(top)
+                .copied()
+                .ok_or(QRError::IndexOutOfBounds(top, forest_len))?;
+            let top_label = forest
+                .vertex_label(top_index)?
+                .ok_or(QRError::NodeNoLabel(top_index))?;
+
+            if top_label.is_packed() {
+                if forest.is_empty_node(top_index)? {
+                    if self.packs.is_empty() {
+                        return Ok(State::Finish);
+                    } else {
+                        return Ok(State::ChangeBranch);
+                    }
+                }
+
+                top = forest.children_of(top_index)?.next().unwrap();
+
+                continue 'search;
+            }
+
+            let query_label = self
+                .query_dfs
+                .get(qtop)
+                .copied()
+                .and_then(|ind| self.query.label_of_vertex(ind).ok());
+
+            if query_label
+                .map(|label| label.matches(top_label.label(), top_label.status()))
+                .unwrap_or(true)
+            {
+                if let Some(captured_value) = candidate_captures.get_mut(qtop) {
+                    *captured_value = top;
+                }
+
+                top += 1;
+                qtop += 1;
+
+                if query_empty_p || qtop >= query_len {
+                    self.current.replace(candidate_captures);
+
+                    return Ok(State::Continue);
+                }
+            } else if qtop == 0 {
+                if self.packs.is_empty() {
+                    return Ok(State::Finish);
+                } else {
+                    return Ok(State::ChangeBranch);
+                }
+            } else {
+                qtop -= 1;
+            }
+        }
+    }
+
+    #[inline]
+    fn clear(&mut self) {
+        self.current.take();
+        self.packs.clear();
+    }
+}
+
+#[cfg(test)]
+mod test_forest_query {
+    use std::io::Write;
+
+    use super::*;
+
+    #[test]
+    fn test_read_forest() -> Result<(), Box<dyn std::error::Error>> {
+        let mut file = std::fs::OpenOptions::new().read(true).open("forest.rep")?;
+
+        use std::io::Read;
+
+        let mut source: Vec<u8> = Vec::new();
+
+        file.read_to_end(&mut source)?;
+
+        use crate::bytes::BytesForest;
+
+        let bytes_forest = BytesForest::try_from(source)?;
+
+        dbg!(
+            bytes_forest.nodes_len(),
+            bytes_forest.children_of(0)?.collect::<Vec<_>>()
+        );
+
+        #[cfg(feature = "test-print-viz")]
+        bytes_forest.print_viz("test forest.gv")?;
+
+        let preprocessed = preprocess_forest(bytes_forest, 0)?;
+
+        // dbg!(&preprocessed);
+
+        use std::io::BufRead;
+
+        println!("Input a query:");
+        std::io::stdout().flush()?;
+
+        for line in std::io::stdin().lock().lines() {
+            let line = line?;
+
+            if line == "q" {
+                break;
+            }
+
+            let qtype = if line.starts_with('s') {
+                1u8
+            } else if line.starts_with('n') {
+                2u8
+            } else {
+                0u8
+            };
+
+            let type_str: &str;
+
+            let q = if line.starts_with('s') || line.starts_with('n') {
+                if let Ok(n) = line[1..].parse::<usize>() {
+                    n
+                } else {
+                    continue;
+                }
+            } else if let Ok(n) = line.parse::<usize>() {
+                n
+            } else {
+                continue;
+            };
+
+            let ans = match qtype {
+                0 => {
+                    type_str = "dfs";
+                    format!("{:?}", preprocessed.dfs_next.get(q).copied())
+                }
+                1 => {
+                    type_str = "subtree range";
+                    format!("{:?}", preprocessed.subtree_ranges.get(q))
+                }
+                _ => {
+                    type_str = "next branch";
+                    format!("{:?}", "haha")
+                }
+            };
+
+            println!("answer for {type_str} = {ans}");
+
+            println!("Input a query:");
+            std::io::stdout().flush()?;
+        }
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_preprocess() -> Result<(), Box<dyn std::error::Error>> {
+        Ok(())
+    }
+}
diff --git a/forest/src/iter/mod.rs b/forest/src/iter/mod.rs
index 49aa4f1..5c55943 100644
--- a/forest/src/iter/mod.rs
+++ b/forest/src/iter/mod.rs
@@ -1,6 +1,6 @@
 //! This module defines the expected behaviours of "iterators of
 //! forests", or "walkers of forests".  Moreover, an implementation
-//! that provides querying capabilities is provided.
+//! with querying capabilities is provided.
 
 use super::*;
 
@@ -30,17 +30,62 @@ impl std::fmt::Display for State {
     }
 }
 
+// TODO: We need a trait to define the behaviours of "queryable"
+// forests: the forests with the three following functionalities:
+//
+// 1. Follow nodes in the depth-first order.
+//
+// 2. Detect if a node is packed.
+//
+// 3. Go to the next "world line" in the forest; that is, go up and
+//   stop before the nearest parent packed node and then go to the
+//   next sibling.
+//
+// We shall call this trait "Queryable".
+
+/// This trait defines the behaviours of data types that can be walked
+/// by walkers specified in this module.
+#[allow(dead_code)]
+pub(in crate::iter) trait Queryable {
+    /// Return the index of the currently visited node.
+    ///
+    /// If it is not visiting nodes, then `None` is returned.
+    fn current(&self) -> Option<usize>;
+
+    /// Go to the next node according to the depth-first order.
+    ///
+    /// If and only if no next node is available, return false.
+    fn next(&mut self) -> bool;
+
+    /// Go to the next branch.  If and only if no next branch to
+    /// explore is available, return false.
+    ///
+    /// By going to the next branch we mean to go up and stop just
+    /// before the nearest parent packed node and then go to the next
+    /// sibling, if any.  If no next sibling is available, repeat the
+    /// above process, until an available branch is found, or until
+    /// every branch is exhausted and we return false.
+    fn next_branch(&mut self) -> bool;
+
+    /// Restrict the forest to the subtree rooted under a specified
+    /// node.
+    ///
+    /// Return whether the currently visited node is rooted under the
+    /// specified node or not.
+    fn restrict(&mut self, node: usize) -> bool;
+}
+
 /// An iterator is able to walk any forest, so no forest type is
 /// associated with iterators.  Essentially speaking, an iterator just
 /// provides the abilities to iterate through nodes of forests which
 /// meet certain conditions, in certain orders.
 pub trait ForestIter: Default {
-    /// Return the currenty visited node.
+    /// Return the currenty visited nodes.
     ///
     /// If the iterator is not visiting any node, return `None`.
     ///
     /// The default implementation always returns `None`.
-    fn current(&self) -> Option<usize> {
+    fn current(&self) -> Option<&[usize]> {
         None
     }
 
diff --git a/forest/src/iter/plan.org b/forest/src/iter/plan.org
new file mode 100644
index 0000000..4ef7db1
--- /dev/null
+++ b/forest/src/iter/plan.org
@@ -0,0 +1,187 @@
+#+TITLE: Plan for querying
+#+AUTHOR: बुबोधयिषुः
+#+DATE: <2024-01-11 Jeu 11:40>
+
+* Basic plan
+
+Basically speaking, we are trying to solve a subgraph isomorphism
+problem, but for parse forests only.  Moreover, each pattern node
+specifies some condition that the target node should satisfy, like the
+rule number, the grammar node number et cetera.  This simply acts like
+the filtering algorithms employed by the general subgraph isomorphism
+algorithms.
+
+There is however one difference from the original problem:
+
+#+begin_quote
+Edges might be matched by paths.
+#+end_quote
+
+This means that some edges in the pattern graph could be matched by
+paths in the (directed) target graph.  Moreover, some edges require to
+be siblings of some previously matched edges.  This poses some
+traditionally missing contraints.
+
+But anyways, we can consult the backtracing algorithms that were
+written for the general problem to have an idea how to solve our
+specific variant of the problem.
+
+** Filtering
+
+Our pattern nodes already specify some conditions for the target nodes
+to match.  So it is worth condidering whether or not we employ some
+more powerful filtering techniques, like the iterative strenghthening
+of the subgraph-isomorphism-compatible-labels algorithm.  This can be
+simply adapted to our situation since we have a natural candidate for
+the labels: the pattern nodes have the matching conditions as labels,
+whereas the target nodes have labels that are either grammar nodes or
+rule numbers.  We can define a partial order on this set of labels by
+requiring that a matching condition is less than or equal to a target
+label if and only if that target label satisfies the matching
+condition.
+
+Then we simply apply the iterative strenghthening procedure, by
+forming, from the label \(L\), the new label \(L\times N(L)\), where
+\(N(L)\) is the function that sends each node \(n\) to the /multiset/
+\(\{L(m) \mid \text{Adj}(m, n)\}\), where \(\text{Adj}\) means that
+the two nodes are adjacent.  And we say that a multiset \(m_1\) is less
+than or equal to another multiset \(m_2\) if and only if the bipartite
+graph formed from these two multisets have a matching covering every
+vertex corresponding to elements in \(m_1\).  Note that this bipartite
+graph has as nodes \(N_x^i\) for every \(x\in m_1\), where \(i = 0, \ldots,
+m_1(x)\) and \(m_1(x)\) is the multiplicity of \(x\) in \(m_1\).  Also
+two nodes \(N_x^i\) and \(N_y^j\), where \(x\in m_1\) and \(y\in m_2\) are
+adjacent if and only if \(x\leq y\).  We can prove that if \(L\) is
+subgraph-isomorphism-compatible, then this new label \(L\times N(L)\) is
+also subgraph-isomorphism-compatible.
+
+As a consequence, we can iteratively construct new labels to filter
+out wrong candidate pairs of pattern nodes and target nodes, until
+this procedure does not yield stronger labels, which are guaranteed to
+happen in polynomial time.
+
+** General sketch
+
+Since our starting label is already domain-specific, it might be
+anticipated that this procedure can enhance the performance of the
+algorithm a lot.  This is the general sketch of the algorithm.
+
+* Two difficulties
+
+** When a branch is exhausted, how should we switch to another branch?
+
+To solve this issue, we need a data type with the following methods:
+
+- Follow nodes in the depth-first order
+- Enumerate and switch to a specific "world line" in the forest: that
+  is, a forest is a pack of trees, and we match against each tree
+  separately; when a tree cannot be matched, we switch to the next
+  tree.
+
+All of this should be implemented efficiently, that is, within a
+constant and "small" amount of operations.
+
+For now, the plan is to preprocess the forest to store the information
+about each packed node in advance, and then the world line is defined
+by which child of those packed nodes to use.  Thus to switch world
+line, we just change an array of numbers, which is pretty efficient I
+suppose.
+
+** When a query cannot be matched, how should we fallback to a parent query?
+
+To solve this issue, we need a data type with the following methods:
+
+- Follow nodes in the depth-first order
+- Jump to the parent and restore the querying context that was stored
+  earlier.
+  
+  We can consult the codes for general subgraph isomorphism problems
+  to implement this, since the general structure is the same.
+
+* Steps [1/4]
+
+** TO-THINK Preprocess the target graph
+
+The target needs to be preprocessed so that we know about the world
+lines in advance, and we can generate some supplemental graphs as well
+to facilitate the search.
+
+This may not be required anymore.
+
+** SOLVED Backtrack based on domains
+We need to keep track of the domains.  Can we find better ways to do
+this?
+
+Maybe we can keep track of only the differences of the set of domains,
+where each domain is a bitset.  For the functionality of bit-sets,
+[[https://docs.rs/bitset-core/latest/bitset_core/index.html][this crate]] is quite useful.
+
+Since this is an external crate, we may treat this feature as an
+optionally enabled feature, and by default we use ordinary arrays of
+booleans.  The main advantage, in my opinion, provided by the external
+crate is its SIMD-awareness, so the default implementation should not
+lose exponential performance.
+
+A more serious problem is how to handle the packed nodes, or
+equivalently how to treat different trees packed inside a forest.
+Previously I was thinking of searching each packed tree separately,
+but this is not respecting the shared structures of these trees, so
+that we must waste a lot of efforts to search the same portions
+multiple times.
+
+A better approach could be to ignore all children of packed nodes, and
+treat the grand-chidlren of the packed nodes as their own children.
+This way we do not need to waste efforts, and we need no special
+treatment of different world-lines when searching.
+
+Here is a more formal description of the approach.
+
+- For each pattern vertex, associate a vector of target vertices,
+  initialized as a vector of vertices with the same boolean true
+  label of length the number of nodes of the target graph.
+- Repeat the following steps.
+- Find a pattern vertex with the smallest domain, tiebreaking by
+  randomly choosing one with the higher degree, excluding vertices
+  whose domains are either empty or singleton.
+  + If there is a vertex with empty domain, the pattern graph is not a
+    subgraph of the target graph.
+  + If every vertex has a singleton domain, return the found mapping.
+- Find a "suitable" target vertex with highest degree within the
+  domain of the chosen vertex.
+  + If none exists, then we go back to the previous association of
+    some target vertex to some pattern vertex, and exclude that target
+    vertex from the domain of that pattern vertex, and then repeat the
+    process again.
+    * If no previous association exists, this means we happen to
+      choose one pattern vertex that cannot be mapped to the target
+      graph at the first step, so we can simply stop the algorithm and
+      decide that the pattern graph is not a subgraph of the target
+      graph.
+- Record the mapping from the above pattern to the above target in the
+  list of associations, and filter other domains.
+  + We can filter by the degree.
+  + We can filter by the sequence of degrees of neighbors.
+
+** TODO Rapid restarts and no-good recording
+This is a highly important optional improvement.
+
+** TODO Connect with the rest of the package
+
+
+Now we want to connect the subgraph isomorphism solver with the rest
+of the package.  This requires a mechanism to convert from user-input
+to a pattern graph, and then run the subgraph (isomorphism) solver to
+find a mapping from the pattern graph to the parse-graph.
+
+In addition, there are some other constraints that might be imposed by
+the user, such as to require that certain edges be siblings, or that
+certain edges be non-repeating (either for nodes or for edges) paths
+instead of plain edges.
+
+All these can be implemented at the filtering step as additional
+requirements other than the neighbourhood-preserving constraint.
+
+Also the user may expect the neighbour relation to be directed.  I am
+not sure if this causes any trouble, but I need to deal with this
+correctly to avoid potential disasters.
+
diff --git a/forest/src/iter/query.rs b/forest/src/iter/query.rs
deleted file mode 100644
index aff8eef..0000000
--- a/forest/src/iter/query.rs
+++ /dev/null
@@ -1,372 +0,0 @@
-//! This module implements an iterator of forests, which provides
-//! querying capabilities.
-
-use super::*;
-
-use graph::{error::Error as GError, labelled::SLGraph, Graph};
-
-use graph_macro::Graph;
-
-/// The label of a query graph.
-///
-/// Such a label is responsible for matching forest labels.
-pub trait QueryLabel: GraphLabel {
-    /// The type of labels this matches against.
-    type Obj: GraphLabel;
-
-    /// Determine whether or not a label matches this query label.
-    fn matches(&self, label: Self::Obj, status: ForestLabelType) -> bool {
-        let _ = label;
-        let _ = status;
-
-        false
-    }
-
-    /// Perform the query for one single query label.
-    fn run<F>(
-        &self,
-        forest: &F,
-        node: usize,
-        packs: &mut Vec<(usize, usize)>,
-    ) -> Result<Option<usize>, QRError>
-    where
-        F: LabelGraph<ForestLabel<<Self as QueryLabel>::Obj>>,
-    {
-        let mut dfs_stack: Vec<usize> = vec![node];
-
-        let nodes_len = forest.nodes_len();
-
-        while let Some(node) = dfs_stack.pop() {
-            let label = match forest.vertex_label(node) {
-                Ok(None) => return Err(QRError::NodeNoLabel(node)),
-                Err(_) => return Err(QRError::IndexOutOfBounds(node, nodes_len)),
-                Ok(Some(label)) => label,
-            };
-
-            let status = label.status();
-            let label = label.label();
-
-            if status == ForestLabelType::Packed {
-                packs.push((node, 0));
-
-                dfs_stack.extend(forest.children_of(node).unwrap().take(1));
-
-                continue;
-            }
-
-            if self.matches(label, status) {
-                return Ok(Some(node));
-            }
-
-            dfs_stack.extend(
-                forest
-                    .children_of(node)
-                    .unwrap()
-                    .collect::<Vec<_>>()
-                    .into_iter()
-                    .rev(),
-            );
-        }
-
-        Ok(None)
-    }
-}
-
-/// How should the current query skip nodes
-#[derive(Debug, Clone, Copy, Default, Ord, PartialOrd, Hash, Eq, PartialEq)]
-pub enum SkipState {
-    /// Skip anything in the depth-first search
-    Any,
-    /// Skip siblings
-    Sibling,
-    /// Do not skip
-    #[default]
-    None,
-}
-
-impl std::fmt::Display for SkipState {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        match self {
-            Self::Any => write!(f, "skip anything"),
-            Self::Sibling => write!(f, "skip siblings"),
-            Self::None => write!(f, "no skipping"),
-        }
-    }
-}
-
-impl From<bool> for SkipState {
-    fn from(value: bool) -> Self {
-        if value {
-            Self::Sibling
-        } else {
-            Self::None
-        }
-    }
-}
-
-/// The type that represents queries as dependency graphs of query
-/// labels.
-///
-/// The labels of vertices are query labels that match against forest
-/// labels, whereas the labels of graphs indicate whether or not the
-/// edge should allow skipping siblings before matching.
-#[derive(Debug, Clone, Graph)]
-pub struct Query<T: QueryLabel> {
-    graph: SLGraph<SkipState>,
-    vertex_labels: Vec<T>,
-    root: Option<usize>,
-}
-
-impl<T: QueryLabel> Query<T> {
-    /// Construct a new query from raw fields.
-    pub fn new(graph: SLGraph<SkipState>, vertex_labels: Vec<T>, root: Option<usize>) -> Self {
-        Self {
-            graph,
-            vertex_labels,
-            root,
-        }
-    }
-
-    /// Return the root of this query.
-    ///
-    /// This returns `None` if and only if the graph is empty.
-    pub fn root(&self) -> Option<usize> {
-        self.root
-    }
-
-    /// Return the label of the given node.
-    ///
-    /// This cannot be implemented as a method of the `LabelGraph`
-    /// trait as the type differs from that of the labels.
-    #[inline]
-    pub fn label_of_vertex(&self, node: usize) -> Result<T, GError> {
-        self.vertex_labels
-            .get(node)
-            .copied()
-            .ok_or(GError::IndexOutOfBounds(node, self.vertex_labels.len()))
-    }
-}
-
-impl<T: QueryLabel> Default for Query<T> {
-    fn default() -> Self {
-        let graph = Default::default();
-        let vertex_labels = Vec::new();
-        let root = None;
-
-        Self {
-            graph,
-            vertex_labels,
-            root,
-        }
-    }
-}
-
-// delegation of labelgraph
-
-impl<T: QueryLabel> LabelGraph<SkipState> for Query<T> {
-    type Iter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::Iter<'a>
-    where
-        Self: 'a;
-
-    type TupleIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::TupleIter<'a>
-    where
-        Self: 'a;
-
-    type LabelIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::LabelIter<'a>
-    where
-        Self: 'a;
-
-    type EdgeLabelIter<'a> = <SLGraph<SkipState> as LabelGraph<SkipState>>::EdgeLabelIter<'a>
-    where
-        Self: 'a;
-
-    #[inline]
-    fn edge_label(&self, source: usize, target: usize) -> Result<Self::EdgeLabelIter<'_>, GError> {
-        self.graph.edge_label(source, target)
-    }
-
-    fn find_children_with_label(
-        &self,
-        _node_id: usize,
-        _label: &SkipState,
-    ) -> Result<<Self as LabelGraph<SkipState>>::Iter<'_>, GError> {
-        unimplemented!("Unsupported operation")
-    }
-
-    fn labels_of(&self, _node_id: usize) -> Result<Self::LabelIter<'_>, GError> {
-        unimplemented!("Unsupported operation")
-    }
-
-    #[inline]
-    fn tuples_of(&self, node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
-        self.graph.tuples_of(node_id)
-    }
-
-    #[inline]
-    fn has_edge_label(
-        &self,
-        node_id: usize,
-        label: &SkipState,
-        target: usize,
-    ) -> Result<bool, GError> {
-        self.graph.has_edge_label(node_id, label, target)
-    }
-}
-
-#[derive(Debug, Clone)]
-/// The type of errors that can occur when querying forests.
-pub enum QRError {
-    /// AN index is out of bounds.
-    IndexOutOfBounds(usize, usize),
-    /// A node has no labels while it is required to have one.
-    NodeNoLabel(usize),
-    /// Encounter an invalid error in converting from an error of
-    /// graphs.
-    InvalidGraphError(GError),
-}
-
-impl std::fmt::Display for QRError {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        match self {
-            QRError::IndexOutOfBounds(n, b) => {
-                write!(f, "index {n} out of bound {b}")
-            }
-            QRError::NodeNoLabel(n) => write!(f, "the node {n} has no label but should have one"),
-            QRError::InvalidGraphError(g) => {
-                write!(f, "an invalid error for forest iterators: {g}")
-            }
-        }
-    }
-}
-
-impl std::error::Error for QRError {}
-
-impl From<GError> for QRError {
-    fn from(value: GError) -> Self {
-        match value {
-            GError::IndexOutOfBounds(n, b) => Self::IndexOutOfBounds(n, b),
-            _ => Self::InvalidGraphError(value),
-        }
-    }
-}
-
-/// An implementation of the forest iterator for querying purposes.
-#[derive(Debug, Clone)]
-pub struct QRun<T: QueryLabel> {
-    /// The currently visited node.
-    current: Option<usize>,
-    /// The queries dependency graph.
-    query: Query<T>,
-    /// The packed nodes that can be used to change to another branch
-    /// if needed.
-    packs: Vec<usize>,
-}
-
-impl<T: QueryLabel> Default for QRun<T> {
-    fn default() -> Self {
-        let current = Default::default();
-        let query = Default::default();
-        let packs = Default::default();
-
-        Self {
-            current,
-            query,
-            packs,
-        }
-    }
-}
-
-impl<T: QueryLabel> ForestIter for QRun<T> {
-    #[inline]
-    fn current(&self) -> Option<usize> {
-        self.current
-    }
-
-    type ForestLabel = <T as QueryLabel>::Obj;
-
-    type Item<'a> = Option<&'a [usize]>
-    where
-        Self: 'a;
-
-    type Error = QRError;
-
-    #[allow(unused)]
-    #[allow(clippy::collapsible_if)]
-    fn next<F>(&mut self, forest: &F, root: usize) -> Result<State, Self::Error>
-    where
-        F: LabelGraph<ForestLabel<Self::ForestLabel>>,
-    {
-        let query_root = if let Some(r) = self.query.root() {
-            r
-        } else {
-            return Ok(State::Finish);
-        };
-
-        #[derive(Debug, Clone, Copy)]
-        enum Item<T: Copy + Clone + std::fmt::Debug> {
-            Seen(T),
-            Unseen(T),
-        }
-
-        use Item::*;
-
-        impl<T: Copy + Clone + std::fmt::Debug> Item<T> {
-            fn is_seen(self) -> bool {
-                matches!(self, Self::Seen(_))
-            }
-
-            fn value(&self) -> T {
-                match self {
-                    Seen(v) => *v,
-                    Unseen(v) => *v,
-                }
-            }
-        }
-
-        let query_len = self.query.nodes_len();
-
-        let mut stack: Vec<Item<(SkipState, usize, usize)>> = Vec::with_capacity(query_len);
-
-        stack.push(Unseen((SkipState::Any, root, query_root)));
-
-        let mut parents_stack: Vec<usize> = Vec::with_capacity(query_len);
-
-        let mut match_result = false;
-
-        while let Some(item) = stack.pop() {
-            if item.is_seen() {
-                if match_result {
-                    continue;
-                }
-
-                // fallback
-            }
-
-            let (state, top, qtop) = item.value();
-
-            let top_label = forest.vertex_label(top)?.ok_or(QRError::NodeNoLabel(top))?;
-
-            let query_label = self.query.label_of_vertex(qtop)?;
-
-            if query_label.matches(top_label.label(), top_label.status()) {
-                // Success branch
-
-                match_result = true;
-
-                // TODO: Deal with captures later
-            } else {
-                // Failure branch
-
-                match_result = false;
-            }
-        }
-
-        todo!()
-    }
-
-    #[inline]
-    fn clear(&mut self) {
-        self.current.take();
-        self.packs.clear();
-    }
-}
diff --git a/forest/src/iter/query/matching.rs b/forest/src/iter/query/matching.rs
new file mode 100644
index 0000000..7ff47bb
--- /dev/null
+++ b/forest/src/iter/query/matching.rs
@@ -0,0 +1,338 @@
+//! This file defines some traits and basic implementations for
+//! matching pattern vertices against target vertices.
+
+use crate::{bytes::BytesForest, ForestLabel};
+
+#[allow(unused_imports)]
+use graph::{
+    error::Error as GError, ELabelGraph, Graph, GraphLabel, LabelGraph, PLGraph, SLGBuilder,
+    SLGraph, VLabelGraph,
+};
+
+use graph_macro::Graph;
+
+use grammar::GrammarLabel;
+
+/// An aggregate struct that holds information about a parent as well
+/// as its vertex label and the corresponding edge label.
+#[derive(Debug, Clone, Copy)]
+pub struct EVLabelParent<T: GraphLabel> {
+    parent: usize,
+    // edge_index: usize,
+    parent_label: T,
+    edge_label: EdgeType,
+}
+
+impl<T: GraphLabel> EVLabelParent<T> {
+    /// Construct a new instance.
+    pub fn new(parent: usize, parent_label: T, edge_label: EdgeType) -> Self {
+        Self {
+            parent,
+            parent_label,
+            edge_label,
+        }
+    }
+
+    /// Retrieve the `parent` field.
+    pub fn parent(&self) -> usize {
+        self.parent
+    }
+
+    // /// Retrieve the `edge_index` field.
+    // pub fn edge_index(&self) -> usize {
+    //     self.edge_index
+    // }
+
+    /// Retrieve the `parent_label` field.
+    pub fn parent_label(&self) -> T {
+        self.parent_label
+    }
+
+    /// Retrieve the `edge_label` field.
+    pub fn edge_label(&self) -> EdgeType {
+        self.edge_label
+    }
+}
+
+/// A trait for pattern graphs.
+pub trait Pattern<T: GraphLabel>: VLabelGraph<T> + ELabelGraph<EdgeType> {
+    /// An iterator delegation for parents and vertex and edge labels.
+    type ParentIter<'a>: Iterator<Item = EVLabelParent<T>> + 'a
+    where
+        Self: 'a,
+        T: 'a;
+
+    /// Return an iterator that goes through the parents and edge
+    /// labels.
+    fn parents(&self, node: usize) -> Result<Self::ParentIter<'_>, GError>;
+}
+
+/// A trait for target graphs.
+pub trait Target<T: GraphLabel>: VLabelGraph<T> {
+    /// Return an array of node indices whose labels match the given
+    /// label.
+    fn filter_label(&self, label: T) -> Box<[usize]> {
+        let _label = label;
+
+        (0..self.nodes_len()).collect::<Vec<_>>().into_boxed_slice()
+    }
+}
+
+/// A custom pattern type with a fixed edge label and a generic vertex
+/// label.
+///
+/// This type contains two graphs, with the extra one containing
+/// parentship information.
+#[derive(Debug, Clone, Graph)]
+pub struct SimplePattern<T: GraphLabel> {
+    graph: SLGBuilder<T, EdgeType>,
+    parents_graph: SLGBuilder<T, EdgeType>,
+}
+
+impl<T: GraphLabel> Default for SimplePattern<T> {
+    fn default() -> Self {
+        let graph = SLGBuilder::default();
+        let parents_graph = SLGBuilder::default();
+        Self {
+            graph,
+            parents_graph,
+        }
+    }
+}
+
+impl<T: GraphLabel> SimplePattern<T> {
+    /// Return a new pattern from an underlying graph builder.
+    pub fn new(graph: SLGBuilder<T, EdgeType>) -> Result<Self, GError> {
+        let nodes_len = graph.nodes_len();
+
+        let mut parents_graph: SLGBuilder<T, EdgeType> = SLGBuilder::with_capacity(nodes_len);
+
+        for node in 0..nodes_len {
+            // Nodes must have labels.
+            let label = VLabelGraph::vertex_label(&graph, node)?.unwrap();
+
+            parents_graph.add_vertex(label);
+        }
+
+        for node in 0..nodes_len {
+            for (label, child) in ELabelGraph::edlabren_of(&graph, node)? {
+                parents_graph.add_edge(node, child, label)?;
+            }
+        }
+
+        Ok(Self {
+            graph,
+            parents_graph,
+        })
+    }
+}
+
+impl<T: GraphLabel> VLabelGraph<T> for SimplePattern<T> {
+    type TupleIter<'a> = <SLGBuilder::<T, EdgeType> as VLabelGraph<T>>::TupleIter<'a>
+    where
+        Self: 'a,
+        T: 'a;
+
+    fn verlabren_of(&self, node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
+        self.graph.verlabren_of(node_id)
+    }
+}
+
+impl<T: GraphLabel> LabelGraph<EdgeType> for SimplePattern<T> {
+    type Iter<'a> = std::iter::Empty<usize>
+    where
+        Self: 'a;
+
+    type TupleIter<'a> = <SLGBuilder::<T, EdgeType> as LabelGraph<EdgeType>>::TupleIter<'a>
+    where
+        Self: 'a,
+        T: 'a;
+
+    type LabelIter<'a>  = std::iter::Empty<(&'a EdgeType, std::iter::Empty<usize>)>
+    where
+        Self: 'a,
+        T: 'a;
+
+    type EdgeLabelIter<'a> = std::option::IntoIter<EdgeType>
+    where
+        Self: 'a,
+        T: 'a;
+
+    fn edge_label(&self, source: usize, target: usize) -> Result<Self::EdgeLabelIter<'_>, GError> {
+        LabelGraph::edge_label(&self.graph, source, target)
+    }
+
+    fn find_children_with_label(
+        &self,
+        node_id: usize,
+        label: &EdgeType,
+    ) -> Result<<Self as LabelGraph<EdgeType>>::Iter<'_>, GError> {
+        self.graph.find_children_with_label(node_id, label)
+    }
+
+    fn labels_of(&self, _node_id: usize) -> Result<Self::LabelIter<'_>, GError> {
+        unimplemented!("Unsupported operation!")
+    }
+
+    fn tuples_of(&self, node_id: usize) -> Result<Self::TupleIter<'_>, GError> {
+        self.graph.tuples_of(node_id)
+    }
+
+    fn has_edge_label(
+        &self,
+        node_id: usize,
+        label: &EdgeType,
+        target: usize,
+    ) -> Result<bool, GError> {
+        self.graph.has_edge_label(node_id, label, target)
+    }
+}
+
+/// Iterator delegator for parents.
+#[derive(Debug, Clone)]
+pub struct PatternParents<'a, T>
+where
+    T: GraphLabel + 'a,
+{
+    vertex_iter: <SLGBuilder<T, EdgeType> as VLabelGraph<T>>::TupleIter<'a>,
+    edge_iter: <SLGBuilder<T, EdgeType> as ELabelGraph<EdgeType>>::TupleIter<'a>,
+}
+
+impl<'a, T> PatternParents<'a, T>
+where
+    T: GraphLabel + 'a,
+{
+    /// Construct a new iterator.
+    pub fn new(
+        vertex_iter: <SLGBuilder<T, EdgeType> as VLabelGraph<T>>::TupleIter<'a>,
+        edge_iter: <SLGBuilder<T, EdgeType> as ELabelGraph<EdgeType>>::TupleIter<'a>,
+    ) -> Self {
+        Self {
+            vertex_iter,
+            edge_iter,
+        }
+    }
+}
+
+impl<'a, T: 'a + GraphLabel> Iterator for PatternParents<'a, T> {
+    type Item = EVLabelParent<T>;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        self.vertex_iter.next().and_then(|(vlabel, target)| {
+            self.edge_iter
+                .next()
+                .map(|(elabel, _)| EVLabelParent::new(target, vlabel, elabel))
+        })
+    }
+
+    #[inline]
+    fn size_hint(&self) -> (usize, Option<usize>) {
+        self.vertex_iter.size_hint()
+    }
+}
+
+impl<'a, T: 'a + GraphLabel> ExactSizeIterator for PatternParents<'a, T> {
+    #[inline]
+    fn len(&self) -> usize {
+        self.vertex_iter.len()
+    }
+}
+
+impl<T: GraphLabel> Pattern<T> for SimplePattern<T> {
+    type ParentIter<'a> = PatternParents<'a, T> where T:'a;
+
+    fn parents(&self, node: usize) -> Result<Self::ParentIter<'_>, GError> {
+        let vertex_iter = self.parents_graph.verlabren_of(node)?;
+        let edge_iter = self.parents_graph.edlabren_of(node)?;
+
+        Ok(PatternParents::new(vertex_iter, edge_iter))
+    }
+}
+
+// Implement Target for three types.  Another type "DefaultForest"
+// should implement this trait on its own when the need arises.
+//
+// Note that SimplePattern also implements the Target trait.
+
+impl<T: GraphLabel> Target<T> for SimplePattern<T> {
+    fn filter_label(&self, label: T) -> Box<[usize]> {
+        (0..self.nodes_len())
+            .filter(|node| VLabelGraph::vertex_label(self, *node) == Ok(Some(label)))
+            .collect::<Vec<_>>()
+            .into_boxed_slice()
+    }
+}
+
+impl<T: GraphLabel> Target<T> for PLGraph<T> {
+    fn filter_label(&self, label: T) -> Box<[usize]> {
+        (0..self.nodes_len())
+            .filter(|node| VLabelGraph::vertex_label(self, *node) == Ok(Some(label)))
+            .collect::<Vec<_>>()
+            .into_boxed_slice()
+    }
+}
+
+impl Target<ForestLabel<GrammarLabel>> for BytesForest {
+    fn filter_label(&self, label: ForestLabel<GrammarLabel>) -> Box<[usize]> {
+        (0..self.nodes_len())
+            .filter(|node| VLabelGraph::vertex_label(self, *node) == Ok(Some(label)))
+            .collect::<Vec<_>>()
+            .into_boxed_slice()
+    }
+}
+
+// Edge-filtering types
+
+/// The type of edges that this edge can match against.
+///
+/// For the time being, we can only match either a direct child
+/// (i.e. no extra condition) or a sibling of another node.
+///
+/// In the future, we shall be able to use supplemental graphs to
+/// match edges.
+///
+/// This can be used as [`Graphlabel`][graph::GraphLabel].
+#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Hash, Ord, PartialOrd)]
+pub enum EdgeType {
+    /// Matches a direct child.
+    #[default]
+    Direct,
+    /// Match a sibling of another node.
+    Sibling(usize),
+}
+
+impl std::fmt::Display for EdgeType {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            EdgeType::Direct => write!(f, "direct edge"),
+            EdgeType::Sibling(node) => write!(f, "sibling of {node}"),
+        }
+    }
+}
+
+impl EdgeType {
+    /// Return a Sibling-type edge.
+    pub fn sibling(node: usize) -> Self {
+        Self::Sibling(node)
+    }
+
+    /// Return whether this edge matches.
+    pub fn matches<G: Graph>(&self, g: &G, from: usize, to: usize) -> Result<bool, GError> {
+        match self {
+            EdgeType::Direct => g.has_edge(from, to),
+            EdgeType::Sibling(sibling) => {
+                let mut iter = g.children_of(from)?;
+
+                while let Some(child) = iter.next() {
+                    if child == *sibling {
+                        return Ok(iter.next() == Some(to));
+                    }
+                }
+
+                // sibling not found
+
+                Ok(false)
+            }
+        }
+    }
+}
diff --git a/forest/src/iter/query/mod.rs b/forest/src/iter/query/mod.rs
new file mode 100644
index 0000000..8820704
--- /dev/null
+++ b/forest/src/iter/query/mod.rs
@@ -0,0 +1,643 @@
+//! This module implements an iterator of forests, which provides
+//! querying capabilities.
+
+// use super::*;
+
+#[allow(unused_imports)]
+use graph::{error::Error as GError, labelled::SLGraph, Graph, GraphLabel, VLabelGraph};
+
+use std::collections::HashSet as Set;
+
+pub mod matching;
+
+use matching::{Pattern, Target};
+
+// REVIEW: Find better approaches for our specific graphs, which are
+// almost like trees.
+
+/// Most of my graphs are directed graphs, so a function is used to
+/// simulate undirected graphs, in case needed.
+#[allow(dead_code)]
+fn has_edge_undirected<G: Graph>(
+    graph: &G,
+    from: usize,
+    to: usize,
+    simulationp: bool,
+) -> Result<bool, GError> {
+    if !simulationp {
+        graph.has_edge(from, to)
+    } else {
+        Ok(graph.has_edge(from, to)? || graph.has_edge(to, from)?)
+    }
+}
+
+/// Initialize a vector of true bits, of the given length.
+fn initialize_bits(length: usize) -> Box<[bool]> {
+    std::iter::repeat(true)
+        .take(length)
+        .collect::<Vec<bool>>()
+        .into_boxed_slice()
+}
+
+/// Return an index whose corresponding slice has the smallest length
+/// that is not zero nor one.
+///
+/// If some slice is empty return error.
+///
+/// If every slice has length one, return `Ok(None)`.
+fn choose_pattern(domains: &[Box<[bool]>], candidates: &[bool]) -> Result<Option<usize>, usize> {
+    let mut cands: Vec<usize> = Vec::with_capacity(domains.len());
+
+    let mut current_min = usize::MAX;
+
+    for ((index, _), domain_size) in candidates
+        .iter()
+        .enumerate()
+        .filter(|(_, considerp)| **considerp)
+        .zip(
+            domains
+                .iter()
+                .map(|domain| domain.iter().filter(|element| **element).count()),
+        )
+    {
+        match domain_size {
+            0 => return Err(index),
+            _ => {
+                // dbg!(index, domain_size, current_min);
+                match domain_size.cmp(&current_min) {
+                    std::cmp::Ordering::Less => {
+                        current_min = domain_size;
+
+                        cands.clear();
+                        cands.push(index);
+                    }
+                    std::cmp::Ordering::Equal => {
+                        cands.push(index);
+                    }
+                    std::cmp::Ordering::Greater => (),
+                };
+            }
+        }
+    }
+
+    // TODO: Add random-ness later
+    Ok(cands.first().copied())
+
+    // if let Some(cand) = cands.first() {
+    //
+
+    //     Ok(Some(*cand))
+    // } else {
+    //     Ok(None)
+    // }
+}
+
+/// Find a target vertex within the chosen domain.
+///
+/// PATTERN_INDEX is used to obtain the degree of the pattern vertex,
+/// so as to filter inappropriate target vertices.
+///
+/// MAPPING is used to know whether some neighbour of the pattern
+/// vertex has already been mapped, which information can be used to
+/// filter the target vertices.
+///
+/// Return the first found suitable target vertex with the highest
+/// degree.
+///
+/// If no suitable target vertex is found, return Ok(None).
+fn choose_target<L, P, T>(
+    pattern: &P,
+    target: &T,
+    pattern_index: usize,
+    domain: &[bool],
+    mapping: &[Option<usize>],
+) -> Result<Option<usize>, GError>
+where
+    L: GraphLabel,
+    P: Pattern<L>,
+    T: Target<L>,
+{
+    let mut cands: Vec<usize> = Vec::with_capacity(target.nodes_len());
+
+    let mut current_max_degree = 0usize;
+
+    'target_loop: for target_index in domain
+        .iter()
+        .enumerate()
+        .filter(|(_, v)| **v)
+        .map(|(index, _)| index)
+    {
+        // TODO: Advanced and custom filtering should be added here.
+        //
+        // These are edge-related filtering, so perhaps we can use
+        // edge labels.
+        //
+        // Some filtering types to add:
+        //
+        // 1. Path with a length restriction: (>= or = or <=) a
+        //    threshold.
+        //
+        //    The default condition is to restrict to length-one
+        //    paths.
+        //
+        //    We can preprocess the target graph to associate a fixed
+        //    set of "near" vertices to any vertex, so that we can
+        //    filter within a constant amount of time.
+        //
+        //    Or we can first restrict this type so that we can only
+        //    use small lengths.
+        //
+        // 2. Being a sibling with another edge.
+
+        for (label, elabel, pattern_neighbour) in pattern
+            .verlabren_of(pattern_index)?
+            .zip(pattern.edlabren_of(pattern_index)?)
+            .map(|((label, index), (elabel, _))| (label, elabel, index))
+        {
+            match elabel {
+                matching::EdgeType::Direct => {
+                    if let Some(Some(mapped)) = mapping.get(pattern_neighbour) {
+                        if !target.has_edge(target_index, *mapped)? {
+                            continue 'target_loop;
+                        }
+
+                        let mapped_label = VLabelGraph::vertex_label(target, *mapped)?.unwrap();
+
+                        if mapped_label != label {
+                            continue 'target_loop;
+                        }
+                    }
+                }
+                matching::EdgeType::Sibling(_) => todo!(),
+            }
+        }
+
+        for parent in pattern.parents(pattern_index)? {
+            let label = parent.parent_label();
+            let elabel = parent.edge_label();
+            let parent = parent.parent();
+
+            match elabel {
+                matching::EdgeType::Direct => {
+                    if let Some(Some(mapped)) = mapping.get(parent) {
+                        if !target.has_edge(*mapped, target_index)? {
+                            continue 'target_loop;
+                        }
+
+                        let mapped_label = VLabelGraph::vertex_label(target, *mapped)?.unwrap();
+
+                        if mapped_label != label {
+                            continue 'target_loop;
+                        }
+                    }
+                }
+                _ => todo!(),
+            }
+        }
+
+        // a currently valid target is found
+
+        let target_degree = target.degree(target_index)?;
+
+        if target_degree > current_max_degree {
+            current_max_degree = target_degree;
+
+            cands.clear();
+            cands.push(target_index);
+        } else if target_index == current_max_degree {
+            cands.push(target_index);
+        }
+    }
+
+    // TODO: Add random-ness later.
+    Ok(cands.first().copied())
+    // if let Some(cand) = cands.first() {
+    //     Ok(Some(*cand))
+    // } else {
+    //     Ok(None)
+    // }
+}
+
+/// Detect if the domains contain a hidden hall set.
+///
+/// This is an "intermediate level detection", so maybe some hall sets
+/// are not found.
+fn find_hall_sets(domains: &mut [Box<[bool]>]) -> Option<Box<[usize]>> {
+    let len = domains.len();
+
+    let mut result: Vec<usize> = Vec::with_capacity(len);
+
+    let mut hall_set: Set<usize> = Set::with_capacity(len);
+
+    let mut accumulated_set: Set<usize> = Set::with_capacity(len);
+
+    let mut current_len = 0usize;
+
+    let mut sorted_indices_by_len = (0..len).collect::<Vec<_>>();
+
+    sorted_indices_by_len.sort_by_cached_key(|element| {
+        domains
+            .get(*element)
+            .unwrap()
+            .iter()
+            .filter(|element| **element)
+            .count()
+    });
+
+    for index in sorted_indices_by_len {
+        result.push(index);
+
+        let domain = domains.get_mut(index).unwrap();
+
+        for hall in hall_set.iter() {
+            *domain.get_mut(*hall).unwrap() = false;
+        }
+
+        accumulated_set.extend(
+            domain
+                .iter()
+                .enumerate()
+                .filter(|(_, v)| **v)
+                .map(|(id, _)| id),
+        );
+
+        current_len += 1;
+
+        if !domain.iter().any(|element| *element) || accumulated_set.len() < current_len {
+            return Some(result.into_boxed_slice());
+        }
+
+        if accumulated_set.len() == current_len {
+            hall_set.extend(accumulated_set.iter().copied());
+            accumulated_set.clear();
+            current_len = 0;
+        }
+    }
+
+    None
+}
+
+/// Determine if a pattern graph is a subgraph of a target graph.
+///
+/// If true, then return Ok(Some(array)), where array is an array of
+/// indices, whose n-th element is the index of a target vertex that
+/// is associated to the n-th vertex in the pattern graph.
+///
+/// If false, return Ok(None).
+pub fn is_subgraph<L, P, T>(pattern: P, target: T) -> Result<Option<Box<[usize]>>, GError>
+where
+    L: GraphLabel,
+    P: Pattern<L>,
+    T: Target<L>,
+{
+    let pattern_nodes_len = pattern.nodes_len();
+    let target_nodes_len = target.nodes_len();
+
+    if pattern_nodes_len > target_nodes_len {
+        return Ok(None);
+    }
+
+    let mut candidates: Box<[bool]> = std::iter::repeat(true)
+        .take(pattern_nodes_len)
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
+
+    // Exclude isolated pattern vertices.
+
+    // for (index, cand) in candidates.iter_mut().enumerate() {
+    //     if pattern.is_empty_node(index)? {
+    //         *cand = false;
+    //     }
+    // }
+
+    let mut mapping: Box<[Option<usize>]> = std::iter::repeat(None)
+        .take(pattern_nodes_len)
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
+
+    /// A backtrack-vector: from,  to,    domains-vector,     unvisited vertices
+    type Association = Vec<(usize, usize, Box<[Box<[bool]>]>, Box<[bool]>)>;
+
+    let mut association_vec: Association = Vec::with_capacity(pattern_nodes_len);
+
+    let mut pattern_domains: Box<[Box<[bool]>]> =
+        std::iter::repeat_with(|| initialize_bits(target_nodes_len))
+            .take(pattern_nodes_len)
+            .collect::<Vec<_>>()
+            .into_boxed_slice();
+
+    // TODO: Perhaps consider some supplemental graphs to facilitate
+    // the filtering: for each supplemental graph S, if f is an
+    // injective mapping from P to T, then f induces an injective
+    // mapping from Hom(S, P) to Hom(S, T).  Essentially this is using
+    // the paths information to filter targets.
+
+    // TODO: Some target vertices might be excluded from every domain,
+    // so we can use the degree-filtering to exclude more domains.
+    // This recursive filtering might be useful.
+
+    for (pattern_index, domain) in pattern_domains.iter_mut().enumerate() {
+        let pattern_degree = pattern.degree(pattern_index)?;
+
+        if let Some(label) = pattern.vertex_label(pattern_index)? {
+            for cand in domain.iter_mut() {
+                *cand = false;
+            }
+
+            for true_cand in target.filter_label(label).iter() {
+                *domain.get_mut(*true_cand).unwrap() = true;
+            }
+        }
+
+        for (target_index, target_handle) in domain.iter_mut().enumerate().filter(|(_, b)| **b) {
+            let target_degree = target.degree(target_index)?;
+
+            if pattern_degree > target_degree {
+                *target_handle = false;
+            }
+        }
+    }
+
+    /// A macro for backtracking
+    macro_rules! backtrack {
+        () => {
+            match association_vec.pop() {
+                None => {
+                    return Ok(None);
+                }
+                Some((from, to, prev_domains, prev_candidates)) => {
+                    dbg!((from, to));
+                    *mapping.get_mut(from).unwrap() = None;
+
+                    pattern_domains = prev_domains;
+
+                    *pattern_domains.get_mut(from).unwrap().get_mut(to).unwrap() = false;
+
+                    candidates = prev_candidates;
+                }
+            }
+        };
+    }
+
+    // use a dummy counter for the time being
+    'main_loop: for _ in 0..100 {
+        dbg!(&pattern_domains);
+
+        match choose_pattern(&pattern_domains, &candidates) {
+            Err(index) => {
+                dbg!(format!("{index} proves impossibility"));
+
+                // silence the compiler
+                let _index = index;
+
+                backtrack!();
+            }
+            Ok(None) => {
+                let mut final_mapping: Vec<usize> = Vec::with_capacity(pattern_nodes_len);
+                let mut target_used_p = std::iter::repeat(false)
+                    .take(target_nodes_len)
+                    .collect::<Vec<bool>>()
+                    .into_boxed_slice();
+
+                for domain in pattern_domains.iter() {
+                    // each domain should be a singleton list at this point
+                    let domain_element = domain.iter().position(|x| *x).unwrap();
+
+                    let used_p = target_used_p.get_mut(domain_element).unwrap();
+
+                    if *used_p {
+                        dbg!("repetition found");
+
+                        backtrack!();
+
+                        continue 'main_loop;
+                    }
+
+                    *used_p = true;
+
+                    final_mapping.push(domain_element);
+                }
+
+                return Ok(Some(final_mapping.into_boxed_slice()));
+            }
+            Ok(Some(index)) => {
+                dbg!(index);
+
+                let chosen_domain = pattern_domains.get(index).unwrap();
+
+                let choice = choose_target(&pattern, &target, index, chosen_domain, &mapping)?;
+
+                match choice {
+                    None => {
+                        backtrack!();
+                    }
+                    Some(choice) => {
+                        dbg!((index, choice));
+
+                        // record this mapping
+
+                        *mapping.get_mut(index).unwrap() = Some(choice);
+
+                        let cloned_domains = pattern_domains.to_vec().into_boxed_slice();
+
+                        let cloned_candidates = candidates.to_vec().into_boxed_slice();
+
+                        association_vec.push((index, choice, cloned_domains, cloned_candidates));
+
+                        // Do not consider this index again.
+
+                        *candidates.get_mut(index).unwrap() = false;
+
+                        // update domains
+
+                        for domain in pattern_domains.iter_mut() {
+                            *domain.get_mut(choice).unwrap() = false;
+                        }
+
+                        let domain = pattern_domains.get_mut(index).unwrap();
+
+                        for domain_ref in domain.iter_mut() {
+                            *domain_ref = false;
+                        }
+
+                        *domain.get_mut(choice).unwrap() = true;
+
+                        // all-different propagation
+
+                        let hall_detection = find_hall_sets(&mut pattern_domains);
+
+                        if let Some(fail_set) = hall_detection {
+                            dbg!(format!("a hall set is found: {fail_set:?}"));
+
+                            backtrack!();
+                        }
+
+                        // TODO: Optionally add degree filtering?
+                    }
+                }
+            }
+        }
+    }
+
+    dbg!("loop time out!");
+
+    Ok(None)
+}
+
+// TODO: Reconstruct the tests from the new types.
+
+#[cfg(test)]
+mod test_subgraph_solver {
+    use super::{
+        matching::{EdgeType, SimplePattern},
+        *,
+    };
+
+    #[allow(unused_imports)]
+    use graph::{Builder, SLGBuilder};
+
+    //     /// Simulate undirected graphs by adding edges in both directions.
+    //     fn add_edge<B>(
+    //         builder: &mut B,
+    //         from: usize,
+    //         to: usize,
+    //         label: <B>::Label,
+    //     ) -> Result<(), Box<dyn std::error::Error>>
+    //     where
+    //         B: Builder,
+    //         <B>::Label: Clone + Copy,
+    //     {
+    //         builder.add_edge(from, to, label)?;
+    //         builder.add_edge(to, from, label)?;
+    //         Ok(())
+    //     }
+
+    #[allow(unused)]
+    fn build_line(len: usize) -> Result<SimplePattern<usize>, Box<dyn std::error::Error>> {
+        let mut builder = SLGBuilder::<usize, EdgeType>::default();
+
+        for i in 0..len {
+            builder.add_vertex(i);
+        }
+
+        for i in 0..(len - 1) {
+            builder.add_edge(i, i + 1, EdgeType::Direct)?;
+        }
+
+        Ok(SimplePattern::new(builder)?)
+    }
+
+    #[allow(unused)]
+    fn build_star(len: usize) -> Result<SimplePattern<usize>, Box<dyn std::error::Error>> {
+        let mut builder = SLGBuilder::<usize, EdgeType>::default();
+
+        for i in 0..len {
+            builder.add_vertex(i);
+        }
+
+        for i in 1..len {
+            builder.add_edge(i, 0, EdgeType::default())?;
+        }
+
+        Ok(SimplePattern::new(builder)?)
+    }
+
+    #[allow(unused)]
+    fn build_triangle_plus_one_edge() -> Result<SimplePattern<usize>, Box<dyn std::error::Error>> {
+        let mut builder = SLGBuilder::<usize, EdgeType>::default();
+
+        for i in 0..4 {
+            builder.add_vertex(i);
+        }
+
+        for i in 1..=3 {
+            builder.add_edge(i, 0, Default::default())?;
+        }
+
+        builder.add_edge(1, 2, EdgeType::default())?;
+
+        Ok(SimplePattern::new(builder)?)
+    }
+
+    #[allow(unused)]
+    fn build_two_triangles_connected() -> Result<SimplePattern<usize>, Box<dyn std::error::Error>> {
+        let mut builder = SLGBuilder::<usize, EdgeType>::default();
+
+        for i in 0..6 {
+            builder.add_vertex(i);
+        }
+
+        for i in 0..=2 {
+            builder.add_edge(i, i + 1, EdgeType::default())?;
+        }
+
+        builder.add_edge(0, 2, EdgeType::default())?;
+
+        for i in 3..=4 {
+            builder.add_edge(i, i + 1, EdgeType::default())?;
+        }
+
+        builder.add_edge(3, 5, EdgeType::default())?;
+
+        Ok(SimplePattern::new(builder)?)
+    }
+
+    #[test]
+    fn test_lines() -> Result<(), Box<dyn std::error::Error>> {
+        let pattern = build_line(4)?;
+        let target = build_line(5)?;
+
+        let result = is_subgraph(pattern, target);
+
+        println!("result = {result:?}");
+
+        match result {
+            Ok(Some(result)) => {
+                assert_eq!(result.to_vec(), [0, 1, 2, 3].to_vec());
+
+                Ok(())
+            }
+            Ok(None) => Err("Subgraph solver thinks this is not a subgraph?!".into()),
+            Err(e) => Err(format!("Subgraph solver failed with error {e}").into()),
+        }
+    }
+
+    #[test]
+    fn test_line_and_star() -> Result<(), Box<dyn std::error::Error>> {
+        let pattern = build_line(4)?;
+        let target = build_star(5)?;
+
+        let result = is_subgraph(pattern, target);
+
+        println!("result = {result:?}");
+
+        match result {
+            Ok(Some(result)) => Err(format!(
+                "Subgraph solver thinks this is a subgraph?!\nMapping = {result:?}"
+            )
+            .into()),
+            Ok(None) => Ok(()),
+            Err(e) => Err(format!("Subgraph solver failed with error {e}").into()),
+        }
+    }
+
+    //     #[test]
+    //     fn test_paper_page_317_example() -> Result<(), Box<dyn std::error::Error>> {
+    //         let pattern = build_triangle_plus_one_edge()?;
+    //         let target = build_two_triangles_connected()?;
+
+    //         let result = is_subgraph(pattern, target);
+
+    //         println!("result = {result:?}");
+
+    //         match result {
+    //             Ok(Some(result)) => {
+    //                 assert_eq!(result.to_vec(), [2, 0, 1, 3].to_vec());
+
+    //                 Ok(())
+    //             }
+    //             Ok(None) => Err("Subgraph solver thinks this is not a subgraph?!".into()),
+    //             Err(e) => Err(format!("Subgraph solver failed with error {e}").into()),
+    //         }
+    //     }
+}