Support Markov order > 1 in defm_motif_parser (#19)

* Initial plan

* Add support for order 2+ Markov models in defm_motif_parser

- Updated regex pattern to match multiple bracketed groups
- Added logic to handle both 2-group (backwards compatible) and m_order+1 group formulas
- Created comprehensive test file 16b-defm-counts-with-formulas.cpp
- All tests passing (72 assertions in 19 test cases)

Co-authored-by: gvegayon <893619+gvegayon@users.noreply.github.com>

* Final verification - all tests passing

Verified:
- All 19 test cases pass with 72 assertions
- Order 2 formulas like "{0y1} > {y1} > {0y1}" now work correctly
- Backwards compatibility maintained for existing 2-group formulas
- counter_generic() integration verified

Co-authored-by: gvegayon <893619+gvegayon@users.noreply.github.com>

* Restore barry.hpp file

* Add comprehensive error testing for defm_motif_parser

Added 15 error test cases to cover all throw statements:
- Transition with m_order=0
- Column out of range (all modes)
- LHS without time when m_order>1
- Row out of range (all modes)
- Duplicate terms (all modes)
- RHS with non-m_order row
- Explicit row index mismatch
- Wrong number of groups
- Intercept with past event
- Invalid formula syntax

All tests passing (87 assertions in 19 test cases)

Co-authored-by: gvegayon <893619+gvegayon@users.noreply.github.com>

* Address PR review comments

- Updated function documentation to describe both 2-group and m_order+1 group modes
- Clarified comment for 2-group mode behavior
- Removed unused arrow_start and arrow_end variables
- Improved comment clarity for RHS bracket check
- Added safety check for group_idx with explanation about covariate handling
- All tests still passing (87 assertions)

Co-authored-by: gvegayon <893619+gvegayon@users.noreply.github.com>

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: gvegayon <893619+gvegayon@users.noreply.github.com>

Author: Copilot

include/barry/models/defm/formula.hpp

@@ -31,9 +31,24 @@
  * 
  * ## Transition effects
  * 
- * Transition effects can be specified using two sets of curly brackets and
- * an greater-than symbol, i.e., `{...} > {...}`. The first set of brackets,
- * which we call LHS, can only hold `row id` that are less than `m_order`.
+ * Transition effects can be specified using curly brackets separated by
+ * greater-than symbols ('>').
+ * 
+ * **Two-group mode (backwards compatible):** `{...} > {...}`
+ * - First group (LHS): Variables at times 0 to m_order-1. When m_order > 1,
+ *   row indices must be explicitly specified.
+ * - Second group (RHS): Variables at time m_order. Row indices can be omitted
+ *   and will default to m_order.
+ * 
+ * **Multi-group mode (explicit):** `{...} > {...} > ... > {...}` (m_order+1 groups)
+ * - Each group corresponds to a specific time point (0, 1, ..., m_order).
+ * - Row indices can be omitted and will be inferred from group position.
+ * - If specified, row indices must match the group position.
+ * 
+ * Examples:
+ * - Order 1: `{y0_0} > {y0_1}` or `{y0_0} > {y0}` (both valid)
+ * - Order 2: `{y0_0} > {y0}` (2-group, implicit final time)
+ * - Order 2: `{y0_0} > {y0_1} > {y0_2}` (3-group, all explicit)
  * 
  * 
  * @param formula A string specifying the motif formula (see details).
@@ -63,9 +78,10 @@ inline void defm_motif_parser(
         std::string("\\{\\s*[01]?y[0-9]+(_[0-9]+)?(\\s*,\\s*[01]?y[0-9]+(_[0-9]+)?)*\\s*\\}") +
         std::string("(\\s*x\\s*[^\\s]+([(].+[)])?\\s*)?")
         );
+    // Updated pattern to match one or more bracketed groups separated by '>'
     std::regex pattern_transition(
-        std::string("\\{\\s*[01]?y[0-9]+(_[0-9]+)?(\\s*,\\s*[01]?y[0-9]+(_[0-9]+)?)*\\}\\s*(>)\\s*") +
         std::string("\\{\\s*[01]?y[0-9]+(_[0-9]+)?(\\s*,\\s*[01]?y[0-9]+(_[0-9]+)?)*\\s*\\}") +
+        std::string("(\\s*>\\s*\\{\\s*[01]?y[0-9]+(_[0-9]+)?(\\s*,\\s*[01]?y[0-9]+(_[0-9]+)?)*\\s*\\})+") +
         std::string("(\\s*x\\s*[^\\s]+([(].+[)])?\\s*)?")
         );
 
@@ -103,79 +119,187 @@ inline void defm_motif_parser(
 
         }
 
-        // Will indicate where the arrow is located at
-        size_t arrow_position = match.position(4u);
+        // Find all bracketed groups to determine which time point each variable belongs to
+        std::regex bracket_pattern("\\{[^}]+\\}");
+        std::vector<std::pair<size_t, size_t>> bracket_ranges; // start, end positions
+        
+        auto brackets_begin = std::sregex_iterator(formula.begin(), formula.end(), bracket_pattern);
+        for (auto i = brackets_begin; i != empty; ++i)
+            bracket_ranges.push_back({i->position(), i->position() + i->length()});
 
-        // This pattern will match 
-        std::regex pattern("(0?)y([0-9]+)(_([0-9]+))?");
-
-        auto iter = std::sregex_iterator(formula.begin(), formula.end(), pattern);
+        size_t num_groups = bracket_ranges.size();
 
-        for (auto i = iter; i != empty; ++i)
+        // For backwards compatibility, allow 2 groups (original behavior) or m_order+1 groups (new behavior)
+        if (num_groups == 2)
         {
+            // Two-group mode (backwards compatible):
+            // - First group: variables at time 0 to m_order-1 (must have explicit row when m_order > 1)
+            // - Second group: variables at time m_order (row can be implicit or explicit)
 
-            // Baseline position
-            size_t current_location = i->position(0u);
-
-            // First value true/false
-            bool is_positive = true;
-            if (i->operator[](1u).str() == "0")
-                is_positive = false;
+            // This pattern will match 
+            std::regex pattern("(0?)y([0-9]+)(_([0-9]+))?");
 
-            // Variable position
-            size_t y_col = std::stoul(i->operator[](2u).str());
-            if (y_col >= y_ncol)
-                throw std::logic_error("The proposed column is out of range.");
+            auto iter = std::sregex_iterator(formula.begin(), formula.end(), pattern);
 
-            // Time location
-            size_t y_row;
-            std::string tmp_str = i->operator[](4u).str();
-            if (m_order > 1)
+            for (auto i = iter; i != empty; ++i)
             {
-                // If missing, we replace with the location 
-                if (tmp_str == "")
+
+                // Baseline position
+                size_t current_location = i->position(0u);
+
+                // First value true/false
+                bool is_positive = true;
+                if (i->operator[](1u).str() == "0")
+                    is_positive = false;
+
+                // Variable position
+                size_t y_col = std::stoul(i->operator[](2u).str());
+                if (y_col >= y_ncol)
+                    throw std::logic_error("The proposed column is out of range.");
+
+                // Time location
+                size_t y_row;
+                std::string tmp_str = i->operator[](4u).str();
+                if (m_order > 1)
                 {
+                    // If missing, we replace with the location 
+                    if (tmp_str == "")
+                    {
 
-                    if (current_location > arrow_position)
-                        y_row = m_order;
+                        if (current_location >= bracket_ranges[1].first)
+                            y_row = m_order;
+                        else
+                            throw std::logic_error("LHS of transition must specify time when m_order > 1");
+
+                    } else
+                        y_row = std::stoul(tmp_str);
+
+                    if (y_row > m_order)
+                        throw std::logic_error("The proposed row is out of range.");
+
+
+                } else {
+
+                    // If missing, we replace with the location 
+                    if (tmp_str != "")
+                        y_row = std::stoul(tmp_str);
                     else
-                        throw std::logic_error("LHS of transition must specify time when m_order > 1");
+                        y_row = (current_location < bracket_ranges[1].first ? 0u: 1u);
 
-                } else
-                    y_row = std::stoul(tmp_str);
+                }
 
-                if (y_row > m_order)
-                    throw std::logic_error("The proposed row is out of range.");
+                if (selected[y_col * (m_order + 1) + y_row])
+                    throw std::logic_error(
+                        "The term " + i->str() + " shows more than once in the formula.");
+
+                // Only variables at time m_order can be in the RHS (second bracketed group)
+                if ((current_location >= bracket_ranges[1].first) && (y_row != m_order))
+                    throw std::logic_error(
+                        "Only the row " + std::to_string(m_order) +
+                        " can be specified at the RHS of the motif."
+                        );
+
+                selected[y_col * (m_order + 1) + y_row] = true;
+
+                locations.push_back(y_col * (m_order + 1) + y_row);
+                signs.push_back(is_positive);
+                
+
+            }
+        }
+        else if (num_groups == m_order + 1)
+        {
+            // New behavior: each group corresponds to a time point (0 to m_order)
+            
+            // This pattern will match 
+            std::regex pattern("(0?)y([0-9]+)(_([0-9]+))?");
 
+            auto iter = std::sregex_iterator(formula.begin(), formula.end(), pattern);
 
-            } else {
+            for (auto i = iter; i != empty; ++i)
+            {
 
-                // If missing, we replace with the location 
+                // Baseline position
+                size_t current_location = i->position(0u);
+
+                // First value true/false
+                bool is_positive = true;
+                if (i->operator[](1u).str() == "0")
+                    is_positive = false;
+
+                // Variable position
+                size_t y_col = std::stoul(i->operator[](2u).str());
+                if (y_col >= y_ncol)
+                    throw std::logic_error("The proposed column is out of range.");
+
+                // Determine which bracketed group this variable belongs to
+                // Note: The regex pattern ensures all variables are within bracketed groups,
+                // as the pattern only matches content inside brackets. Variables in covariate
+                // expressions (after 'x') are handled separately by pattern_conditional above.
+                size_t group_idx = 0;
+                bool found_group = false;
+                for (size_t g = 0; g < bracket_ranges.size(); ++g)
+                {
+                    if (current_location >= bracket_ranges[g].first && 
+                        current_location < bracket_ranges[g].second)
+                    {
+                        group_idx = g;
+                        found_group = true;
+                        break;
+                    }
+                }
+                
+                // Safety check: ensure the variable was found in a bracketed group
+                // This should never happen given the regex, but verify for safety
+                if (!found_group)
+                    throw std::logic_error(
+                        "Internal error: variable " + i->str() + 
+                        " not found within any bracketed group.");
+
+                // Time location
+                size_t y_row;
+                std::string tmp_str = i->operator[](4u).str();
+                
+                // If row is explicitly specified, use it; otherwise infer from group position
                 if (tmp_str != "")
+                {
                     y_row = std::stoul(tmp_str);
+                    
+                    // Validate that explicit row matches the expected group position
+                    if (y_row != group_idx)
+                        throw std::logic_error(
+                            "Explicit row index " + std::to_string(y_row) + 
+                            " does not match the position in the formula (group " + 
+                            std::to_string(group_idx) + ").");
+                }
                 else
-                    y_row = (current_location < arrow_position ? 0u: 1u);
-
-            }
+                {
+                    // Infer row from bracketed group position
+                    y_row = group_idx;
+                }
 
-            if (selected[y_col * (m_order + 1) + y_row])
-                throw std::logic_error(
-                    "The term " + i->str() + " shows more than once in the formula.");
+                if (y_row > m_order)
+                    throw std::logic_error("The proposed row is out of range.");
 
-            // Only the end of the chain can be located at position after the
-            // arrow
-            if ((current_location > arrow_position) && (y_row != m_order))
-                throw std::logic_error(
-                    "Only the row " + std::to_string(m_order) +
-                    " can be specified at the RHS of the motif."
-                    );
+                if (selected[y_col * (m_order + 1) + y_row])
+                    throw std::logic_error(
+                        "The term " + i->str() + " shows more than once in the formula.");
 
-            selected[y_col * (m_order + 1) + y_row] = true;
+                selected[y_col * (m_order + 1) + y_row] = true;
 
-            locations.push_back(y_col * (m_order + 1) + y_row);
-            signs.push_back(is_positive);
-            
+                locations.push_back(y_col * (m_order + 1) + y_row);
+                signs.push_back(is_positive);
+                
 
+            }
+        }
+        else
+        {
+            throw std::logic_error(
+                "For a Markov model of order " + std::to_string(m_order) + 
+                ", transition formulas must have either 2 bracketed groups " +
+                "(for backwards compatibility) or exactly " + std::to_string(m_order + 1) + 
+                " bracketed groups (found " + std::to_string(num_groups) + ").");
         }
 
         return;