sim_lib_discrete_graph/
alignment.rs

1//! Certified dynamic-time-warp and edit alignment.
2
3use core::cmp::Ordering;
4
5use crate::{
6    AlgorithmControl, AlgorithmInterrupt, AlgorithmReceipt, FiniteCost, GraphError, NeverInterrupt,
7    control::WorkMeter,
8    cost::{add, compare, validate},
9};
10
11mod verify;
12
13pub use verify::verify_alignment;
14
15/// Window limiting which prefix-pair cells an alignment may visit.
16#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
17pub enum AlignmentWindow {
18    /// Evaluate the full prefix grid.
19    #[default]
20    Unbounded,
21    /// Permit cells whose prefix indices differ by at most `radius`.
22    Radius(usize),
23}
24
25impl AlignmentWindow {
26    fn contains(self, left: usize, right: usize) -> bool {
27        match self {
28            Self::Unbounded => true,
29            Self::Radius(radius) => left.abs_diff(right) <= radius,
30        }
31    }
32}
33
34/// Endpoint semantics for sequence alignment.
35#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
36pub enum AlignmentBoundary {
37    /// Align both complete sequences.
38    #[default]
39    Global,
40    /// Align the complete left query to a contiguous region of the right
41    /// sequence. Right-side prefix and suffix costs are free.
42    Subsequence,
43}
44
45/// Whether to retain the full proof table or only the final rolling row.
46#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
47pub enum AlignmentMemory {
48    /// Retain backpointers and return the full alignment path.
49    #[default]
50    Full,
51    /// Retain two score rows while solving and return score-only evidence.
52    RollingScoreOnly,
53}
54
55/// Costs for advancing only one side of an edit alignment.
56#[derive(Clone, Debug, PartialEq, Eq)]
57pub struct GapPolicy<C> {
58    /// Cost of consuming a left item without a right item.
59    pub delete: C,
60    /// Cost of consuming a right item without a left item.
61    pub insert: C,
62}
63
64impl<C> GapPolicy<C> {
65    /// Builds a gap policy with explicit deletion and insertion costs.
66    pub fn new(delete: C, insert: C) -> Self {
67        Self { delete, insert }
68    }
69}
70
71/// Window, gap, boundary, and memory policy for dynamic alignment.
72#[derive(Clone, Debug, PartialEq, Eq)]
73pub struct DtwPolicy<C> {
74    /// Admissible prefix-grid cells.
75    pub window: AlignmentWindow,
76    /// One-sided step costs.
77    pub gaps: GapPolicy<C>,
78    /// Required start/end coverage.
79    pub boundary: AlignmentBoundary,
80    /// Retained proof material.
81    pub memory: AlignmentMemory,
82}
83
84impl<C> DtwPolicy<C> {
85    /// Builds a global, unbounded, full-memory alignment policy.
86    pub fn new(gaps: GapPolicy<C>) -> Self {
87        Self {
88            window: AlignmentWindow::Unbounded,
89            gaps,
90            boundary: AlignmentBoundary::Global,
91            memory: AlignmentMemory::Full,
92        }
93    }
94
95    /// Returns a copy with a different alignment window.
96    pub fn with_window(mut self, window: AlignmentWindow) -> Self {
97        self.window = window;
98        self
99    }
100
101    /// Returns a copy with different endpoint semantics.
102    pub fn with_boundary(mut self, boundary: AlignmentBoundary) -> Self {
103        self.boundary = boundary;
104        self
105    }
106
107    /// Returns a copy with a different memory policy.
108    pub fn with_memory(mut self, memory: AlignmentMemory) -> Self {
109        self.memory = memory;
110        self
111    }
112}
113
114/// Stable predecessor move in a full alignment certificate.
115#[derive(Clone, Copy, Debug, PartialEq, Eq)]
116pub enum AlignmentMove {
117    /// Consume one item from each sequence.
118    Match,
119    /// Consume one left item.
120    Delete,
121    /// Consume one right item.
122    Insert,
123}
124
125/// One reachable prefix-grid cell in a full alignment certificate.
126#[derive(Clone, Debug, PartialEq, Eq)]
127pub struct AlignmentCell<C> {
128    /// Minimum cost of reaching this cell.
129    pub total_cost: C,
130    /// Stable predecessor move, or `None` at a permitted free start.
131    pub predecessor: Option<AlignmentMove>,
132    /// Cost charged by the predecessor move.
133    pub step_cost: Option<C>,
134}
135
136/// Optimality evidence retained under the selected memory policy.
137#[derive(Clone, Debug, PartialEq, Eq)]
138pub enum AlignmentCertificate<C> {
139    /// Full prefix-grid Bellman table with backpointers.
140    Full {
141        /// Row-major `(left prefix, right prefix)` cells.
142        cells: Vec<Vec<Option<AlignmentCell<C>>>>,
143    },
144    /// Final score row reproduced by rolling-memory evaluation.
145    Rolling {
146        /// Cost at every right prefix after consuming the complete left input.
147        final_row: Vec<Option<C>>,
148        /// Stable selected right endpoint.
149        endpoint: usize,
150    },
151}
152
153/// One operation in a full alignment path.
154#[derive(Clone, Debug, PartialEq, Eq)]
155pub enum AlignmentStep<C> {
156    /// Pair one left and one right item.
157    Match {
158        /// Left input index.
159        left: usize,
160        /// Right input index.
161        right: usize,
162        /// Local pair cost.
163        cost: C,
164    },
165    /// Consume an unmatched left item.
166    Delete {
167        /// Left input index.
168        left: usize,
169        /// Gap cost.
170        cost: C,
171    },
172    /// Consume an unmatched right item.
173    Insert {
174        /// Right input index.
175        right: usize,
176        /// Gap cost.
177        cost: C,
178    },
179}
180
181/// Minimum-cost alignment plus proof and deterministic accounting.
182#[derive(Clone, Debug, PartialEq, Eq)]
183pub struct Alignment<C> {
184    /// Exact minimum score.
185    pub score: C,
186    /// Stable full path, absent in rolling score-only mode.
187    pub steps: Option<Vec<AlignmentStep<C>>>,
188    /// Proof material selected by the memory policy.
189    pub certificate: AlignmentCertificate<C>,
190    /// Cell/edge work accounting.
191    pub receipt: AlgorithmReceipt,
192}
193
194/// Aligns two sequences under a window, gap, boundary, and memory policy.
195pub fn dynamic_time_warp<T, C: FiniteCost>(
196    left: &[T],
197    right: &[T],
198    local_cost: impl Fn(&T, &T) -> C,
199    policy: DtwPolicy<C>,
200) -> Result<Alignment<C>, GraphError> {
201    dynamic_time_warp_with_control(
202        left,
203        right,
204        local_cost,
205        policy,
206        &AlgorithmControl::default(),
207        &NeverInterrupt,
208    )
209}
210
211/// Aligns two sequences under explicit work and cancellation control.
212pub fn dynamic_time_warp_with_control<T, C: FiniteCost>(
213    left: &[T],
214    right: &[T],
215    local_cost: impl Fn(&T, &T) -> C,
216    policy: DtwPolicy<C>,
217    control: &AlgorithmControl,
218    interrupt: &dyn AlgorithmInterrupt,
219) -> Result<Alignment<C>, GraphError> {
220    validate_policy(&policy)?;
221    let memory = peak_memory(left.len(), right.len(), policy.memory)?;
222    let mut meter = WorkMeter::new(control, interrupt, memory)?;
223    let computed = match policy.memory {
224        AlignmentMemory::Full => {
225            let (cells, stats) = full_table(left, right, &local_cost, &policy, Some(&mut meter))?;
226            let endpoint = select_endpoint(&cells, policy.boundary)?;
227            let score = cells[left.len()][endpoint]
228                .as_ref()
229                .expect("selected endpoint is reachable")
230                .total_cost
231                .clone();
232            let steps = reconstruct(&cells, left.len(), endpoint, policy.boundary)?;
233            Computed {
234                score,
235                steps: Some(steps),
236                certificate: AlignmentCertificate::Full { cells },
237                stats,
238            }
239        }
240        AlignmentMemory::RollingScoreOnly => {
241            let (final_row, stats) =
242                rolling_row(left, right, &local_cost, &policy, Some(&mut meter))?;
243            let endpoint = select_rolling_endpoint(&final_row, policy.boundary, right.len())?;
244            let score = final_row[endpoint]
245                .as_ref()
246                .expect("selected endpoint is reachable")
247                .clone();
248            Computed {
249                score,
250                steps: None,
251                certificate: AlignmentCertificate::Rolling {
252                    final_row,
253                    endpoint,
254                },
255                stats,
256            }
257        }
258    };
259    let receipt = meter.finish();
260    debug_assert_eq!(receipt.cells, computed.stats.cells);
261    debug_assert_eq!(receipt.edges, computed.stats.edges);
262    Ok(Alignment {
263        score: computed.score,
264        steps: computed.steps,
265        certificate: computed.certificate,
266        receipt,
267    })
268}
269
270struct Computed<C> {
271    score: C,
272    steps: Option<Vec<AlignmentStep<C>>>,
273    certificate: AlignmentCertificate<C>,
274    stats: Stats,
275}
276
277#[derive(Clone, Copy, Debug, Default)]
278struct Stats {
279    cells: u64,
280    edges: u64,
281}
282
283type Table<C> = Vec<Vec<Option<AlignmentCell<C>>>>;
284
285fn full_table<T, C: FiniteCost>(
286    left: &[T],
287    right: &[T],
288    local_cost: &impl Fn(&T, &T) -> C,
289    policy: &DtwPolicy<C>,
290    mut meter: Option<&mut WorkMeter<'_>>,
291) -> Result<(Table<C>, Stats), GraphError> {
292    let rows = left
293        .len()
294        .checked_add(1)
295        .ok_or_else(|| GraphError::WeightOverflow("alignment rows".to_owned()))?;
296    let columns = right
297        .len()
298        .checked_add(1)
299        .ok_or_else(|| GraphError::WeightOverflow("alignment columns".to_owned()))?;
300    let mut cells = vec![vec![None; columns]; rows];
301    let mut stats = Stats::default();
302    for i in 0..rows {
303        for j in 0..columns {
304            if !policy.window.contains(i, j) {
305                continue;
306            }
307            charge_cell(&mut meter, &mut stats)?;
308            cells[i][j] = compute_cell(
309                i,
310                j,
311                left,
312                right,
313                local_cost,
314                policy,
315                |row, column| cells[row][column].clone(),
316                &mut meter,
317                &mut stats,
318            )?;
319        }
320    }
321    Ok((cells, stats))
322}
323
324#[allow(clippy::too_many_arguments)]
325fn compute_cell<C: FiniteCost, T>(
326    i: usize,
327    j: usize,
328    left: &[T],
329    right: &[T],
330    local_cost: &impl Fn(&T, &T) -> C,
331    policy: &DtwPolicy<C>,
332    lookup: impl Fn(usize, usize) -> Option<AlignmentCell<C>>,
333    meter: &mut Option<&mut WorkMeter<'_>>,
334    stats: &mut Stats,
335) -> Result<Option<AlignmentCell<C>>, GraphError> {
336    if i == 0 && (j == 0 || policy.boundary == AlignmentBoundary::Subsequence) {
337        return Ok(Some(AlignmentCell {
338            total_cost: C::zero(),
339            predecessor: None,
340            step_cost: None,
341        }));
342    }
343    let mut best: Option<(C, AlignmentMove, C)> = None;
344    if i > 0 && j > 0 {
345        charge_edge(meter, stats)?;
346        if let Some(previous) = lookup(i - 1, j - 1) {
347            let cost = local_cost(&left[i - 1], &right[j - 1]);
348            validate_non_negative(&cost, "alignment local cost")?;
349            let total = add(&previous.total_cost, &cost, "alignment match")?;
350            choose(&mut best, total, AlignmentMove::Match, cost)?;
351        }
352    }
353    if i > 0 {
354        charge_edge(meter, stats)?;
355        if let Some(previous) = lookup(i - 1, j) {
356            let total = add(
357                &previous.total_cost,
358                &policy.gaps.delete,
359                "alignment deletion",
360            )?;
361            choose(
362                &mut best,
363                total,
364                AlignmentMove::Delete,
365                policy.gaps.delete.clone(),
366            )?;
367        }
368    }
369    if j > 0 {
370        charge_edge(meter, stats)?;
371        if let Some(previous) = lookup(i, j - 1) {
372            let total = add(
373                &previous.total_cost,
374                &policy.gaps.insert,
375                "alignment insertion",
376            )?;
377            choose(
378                &mut best,
379                total,
380                AlignmentMove::Insert,
381                policy.gaps.insert.clone(),
382            )?;
383        }
384    }
385    Ok(
386        best.map(|(total_cost, predecessor, step_cost)| AlignmentCell {
387            total_cost,
388            predecessor: Some(predecessor),
389            step_cost: Some(step_cost),
390        }),
391    )
392}
393
394fn rolling_row<T, C: FiniteCost>(
395    left: &[T],
396    right: &[T],
397    local_cost: &impl Fn(&T, &T) -> C,
398    policy: &DtwPolicy<C>,
399    mut meter: Option<&mut WorkMeter<'_>>,
400) -> Result<(Vec<Option<C>>, Stats), GraphError> {
401    let columns = right
402        .len()
403        .checked_add(1)
404        .ok_or_else(|| GraphError::WeightOverflow("alignment columns".to_owned()))?;
405    let mut previous = vec![None; columns];
406    let mut stats = Stats::default();
407    for i in 0..=left.len() {
408        let mut current = vec![None; columns];
409        for j in 0..columns {
410            if !policy.window.contains(i, j) {
411                continue;
412            }
413            charge_cell(&mut meter, &mut stats)?;
414            if i == 0 && (j == 0 || policy.boundary == AlignmentBoundary::Subsequence) {
415                current[j] = Some(C::zero());
416                continue;
417            }
418            let mut best: Option<C> = None;
419            if i > 0 && j > 0 {
420                charge_edge(&mut meter, &mut stats)?;
421                if let Some(prior) = &previous[j - 1] {
422                    let cost = local_cost(&left[i - 1], &right[j - 1]);
423                    validate_non_negative(&cost, "alignment local cost")?;
424                    choose_score(&mut best, add(prior, &cost, "alignment match")?)?;
425                }
426            }
427            if i > 0 {
428                charge_edge(&mut meter, &mut stats)?;
429                if let Some(prior) = &previous[j] {
430                    choose_score(
431                        &mut best,
432                        add(prior, &policy.gaps.delete, "alignment deletion")?,
433                    )?;
434                }
435            }
436            if j > 0 {
437                charge_edge(&mut meter, &mut stats)?;
438                if let Some(prior) = &current[j - 1] {
439                    choose_score(
440                        &mut best,
441                        add(prior, &policy.gaps.insert, "alignment insertion")?,
442                    )?;
443                }
444            }
445            current[j] = best;
446        }
447        previous = current;
448    }
449    Ok((previous, stats))
450}
451
452fn choose<C: FiniteCost>(
453    best: &mut Option<(C, AlignmentMove, C)>,
454    total: C,
455    movement: AlignmentMove,
456    step_cost: C,
457) -> Result<(), GraphError> {
458    let replace = match best {
459        Some((current, _, _)) => {
460            compare(&total, current, "alignment candidate ordering")? == Ordering::Less
461        }
462        None => true,
463    };
464    if replace {
465        *best = Some((total, movement, step_cost));
466    }
467    Ok(())
468}
469
470fn choose_score<C: FiniteCost>(best: &mut Option<C>, total: C) -> Result<(), GraphError> {
471    let replace = match best {
472        Some(current) => {
473            compare(&total, current, "alignment candidate ordering")? == Ordering::Less
474        }
475        None => true,
476    };
477    if replace {
478        *best = Some(total);
479    }
480    Ok(())
481}
482
483fn select_endpoint<C: FiniteCost>(
484    cells: &Table<C>,
485    boundary: AlignmentBoundary,
486) -> Result<usize, GraphError> {
487    let final_row = cells.last().expect("alignment table has a prefix row");
488    select_cell_endpoint(final_row, boundary)
489}
490
491fn select_cell_endpoint<C: FiniteCost>(
492    row: &[Option<AlignmentCell<C>>],
493    boundary: AlignmentBoundary,
494) -> Result<usize, GraphError> {
495    match boundary {
496        AlignmentBoundary::Global => row
497            .len()
498            .checked_sub(1)
499            .filter(|endpoint| row[*endpoint].is_some())
500            .ok_or(GraphError::Disconnected),
501        AlignmentBoundary::Subsequence => {
502            let mut best: Option<(usize, &C)> = None;
503            for (index, cell) in row.iter().enumerate() {
504                let Some(cell) = cell else {
505                    continue;
506                };
507                let replace = match best {
508                    Some((_, cost)) => {
509                        compare(&cell.total_cost, cost, "alignment endpoint ordering")?
510                            == Ordering::Less
511                    }
512                    None => true,
513                };
514                if replace {
515                    best = Some((index, &cell.total_cost));
516                }
517            }
518            best.map(|(index, _)| index).ok_or(GraphError::Disconnected)
519        }
520    }
521}
522
523fn select_rolling_endpoint<C: FiniteCost>(
524    row: &[Option<C>],
525    boundary: AlignmentBoundary,
526    right_len: usize,
527) -> Result<usize, GraphError> {
528    match boundary {
529        AlignmentBoundary::Global => row
530            .get(right_len)
531            .and_then(Option::as_ref)
532            .map(|_| right_len)
533            .ok_or(GraphError::Disconnected),
534        AlignmentBoundary::Subsequence => {
535            let mut best: Option<(usize, &C)> = None;
536            for (index, cost) in row.iter().enumerate() {
537                let Some(cost) = cost else {
538                    continue;
539                };
540                let replace = match best {
541                    Some((_, current)) => {
542                        compare(cost, current, "alignment endpoint ordering")? == Ordering::Less
543                    }
544                    None => true,
545                };
546                if replace {
547                    best = Some((index, cost));
548                }
549            }
550            best.map(|(index, _)| index).ok_or(GraphError::Disconnected)
551        }
552    }
553}
554
555fn reconstruct<C: FiniteCost>(
556    cells: &Table<C>,
557    mut i: usize,
558    mut j: usize,
559    boundary: AlignmentBoundary,
560) -> Result<Vec<AlignmentStep<C>>, GraphError> {
561    let mut reversed = Vec::new();
562    loop {
563        let cell = cells[i][j].as_ref().ok_or_else(|| {
564            GraphError::CertificateInvalid("alignment path visits an unreachable cell".to_owned())
565        })?;
566        let Some(movement) = cell.predecessor else {
567            if i == 0 && (j == 0 || boundary == AlignmentBoundary::Subsequence) {
568                break;
569            }
570            return Err(GraphError::CertificateInvalid(
571                "alignment path ends at an illegal free boundary".to_owned(),
572            ));
573        };
574        let cost = cell.step_cost.clone().ok_or_else(|| {
575            GraphError::CertificateInvalid("alignment step has no cost".to_owned())
576        })?;
577        match movement {
578            AlignmentMove::Match => {
579                i -= 1;
580                j -= 1;
581                reversed.push(AlignmentStep::Match {
582                    left: i,
583                    right: j,
584                    cost,
585                });
586            }
587            AlignmentMove::Delete => {
588                i -= 1;
589                reversed.push(AlignmentStep::Delete { left: i, cost });
590            }
591            AlignmentMove::Insert => {
592                j -= 1;
593                reversed.push(AlignmentStep::Insert { right: j, cost });
594            }
595        }
596    }
597    reversed.reverse();
598    Ok(reversed)
599}
600
601fn validate_policy<C: FiniteCost>(policy: &DtwPolicy<C>) -> Result<(), GraphError> {
602    validate_non_negative(&policy.gaps.delete, "alignment deletion gap")?;
603    validate_non_negative(&policy.gaps.insert, "alignment insertion gap")
604}
605
606fn validate_non_negative<C: FiniteCost>(cost: &C, context: &str) -> Result<(), GraphError> {
607    validate(cost, context)?;
608    if compare(cost, &C::zero(), context)? == Ordering::Less {
609        return Err(GraphError::Unsupported(format!(
610            "{context} must be non-negative"
611        )));
612    }
613    Ok(())
614}
615
616fn peak_memory(
617    left_len: usize,
618    right_len: usize,
619    memory: AlignmentMemory,
620) -> Result<usize, GraphError> {
621    let columns = right_len
622        .checked_add(1)
623        .ok_or_else(|| GraphError::WeightOverflow("alignment columns".to_owned()))?;
624    match memory {
625        AlignmentMemory::Full => left_len
626            .checked_add(1)
627            .and_then(|rows| rows.checked_mul(columns))
628            .ok_or_else(|| GraphError::WeightOverflow("alignment table cells".to_owned())),
629        AlignmentMemory::RollingScoreOnly => {
630            let rows: usize = if left_len == 0 { 1 } else { 2 };
631            rows.checked_mul(columns)
632                .ok_or_else(|| GraphError::WeightOverflow("alignment rolling rows".to_owned()))
633        }
634    }
635}
636
637fn charge_cell(
638    meter: &mut Option<&mut WorkMeter<'_>>,
639    stats: &mut Stats,
640) -> Result<(), GraphError> {
641    if let Some(meter) = meter.as_deref_mut() {
642        meter.cell()?;
643    }
644    stats.cells = stats
645        .cells
646        .checked_add(1)
647        .ok_or_else(|| GraphError::WeightOverflow("alignment cell count".to_owned()))?;
648    Ok(())
649}
650
651fn charge_edge(
652    meter: &mut Option<&mut WorkMeter<'_>>,
653    stats: &mut Stats,
654) -> Result<(), GraphError> {
655    if let Some(meter) = meter.as_deref_mut() {
656        meter.edge()?;
657    }
658    stats.edges = stats
659        .edges
660        .checked_add(1)
661        .ok_or_else(|| GraphError::WeightOverflow("alignment edge count".to_owned()))?;
662    Ok(())
663}
664
665#[cfg(test)]
666mod tests;
sim_lib_discrete_graph/alignment.rs

sim_lib_discrete_graph/
alignment.rs
