Skip to main content

cosmolkit_core/bio/
ops.rs

1//! Operation contract system for BioStructure.
2//!
3//! Mirrors the `MoleculeOpSpec` / `OpParts` pattern in `ops.rs`, adapted for
4//! the BioStructure hierarchy. Operation bodies must be registered via the
5//! `bio_structure_ops!` macro. `BioOpParts` is wrapper-owned migration and
6//! contract-recording machinery; operation-specific behavior belongs in the
7//! operation body or in bio domain modules, not in `BioOpParts`.
8
9use std::marker::PhantomData;
10
11use crate::{
12    SupportStatus,
13    bio::{AssemblyId, AtomId, BioStructure, BondId, ChainId, EntityId, ModelId, ResidueId},
14};
15
16// ---------------------------------------------------------------------------
17// BioBlockSet — bitmask of mutable blocks
18// ---------------------------------------------------------------------------
19
20#[derive(Debug, Clone, Copy, PartialEq, Eq)]
21pub struct BioBlockSet(u32);
22
23impl BioBlockSet {
24    pub const NONE: Self = Self(0);
25    pub const ATOMS: Self = Self(1 << 0);
26    pub const RESIDUES: Self = Self(1 << 1);
27    pub const CHAINS: Self = Self(1 << 2);
28    pub const ENTITIES: Self = Self(1 << 3);
29    pub const MODELS: Self = Self(1 << 4);
30    pub const COORDINATES: Self = Self(1 << 5);
31    pub const BONDS: Self = Self(1 << 6);
32    pub const ASSEMBLIES: Self = Self(1 << 7);
33    pub const ANNOTATIONS: Self = Self(1 << 8);
34    pub const DERIVED_CACHE: Self = Self(1 << 9);
35    pub const PROPERTIES: Self = Self(1 << 10);
36
37    #[must_use]
38    pub const fn union(self, other: Self) -> Self {
39        Self(self.0 | other.0)
40    }
41
42    #[must_use]
43    pub const fn contains(self, other: Self) -> bool {
44        (self.0 & other.0) == other.0
45    }
46}
47
48// ---------------------------------------------------------------------------
49// BioStateSet — bitmask of structural state that must be handled
50// ---------------------------------------------------------------------------
51
52#[derive(Debug, Clone, Copy, PartialEq, Eq)]
53pub struct BioStateSet(u32);
54
55impl BioStateSet {
56    pub const NONE: Self = Self(0);
57    pub const HIERARCHY: Self = Self(1 << 0);
58    pub const RESIDUE_SPANS: Self = Self(1 << 1);
59    pub const CHAIN_SPANS: Self = Self(1 << 2);
60    pub const MODEL_SPANS: Self = Self(1 << 3);
61    pub const COORDINATE_ALIGNMENT: Self = Self(1 << 4);
62    pub const ENTITY_MAPPING: Self = Self(1 << 5);
63    pub const ALTLOC_GROUPS: Self = Self(1 << 6);
64    pub const ASSEMBLY_REFERENCES: Self = Self(1 << 7);
65    pub const BOND_REFERENCES: Self = Self(1 << 8);
66    pub const SELECTION_PROVENANCE: Self = Self(1 << 9);
67    pub const POLYMER_ANNOTATION: Self = Self(1 << 10);
68    pub const SECONDARY_STRUCTURE: Self = Self(1 << 11);
69
70    #[must_use]
71    pub const fn union(self, other: Self) -> Self {
72        Self(self.0 | other.0)
73    }
74
75    #[must_use]
76    pub const fn contains(self, other: Self) -> bool {
77        (self.0 & other.0) == other.0
78    }
79}
80
81// ---------------------------------------------------------------------------
82// BioDerivedState — bitmask of derived cache entries
83// ---------------------------------------------------------------------------
84
85#[derive(Debug, Clone, Copy, PartialEq, Eq)]
86pub struct BioDerivedState(u64);
87
88impl BioDerivedState {
89    pub const NONE: Self = Self(0);
90    pub const ATOM_INDEX: Self = Self(1 << 0);
91    pub const RESIDUE_INDEX: Self = Self(1 << 1);
92    pub const CHAIN_INDEX: Self = Self(1 << 2);
93    pub const ENTITY_INDEX: Self = Self(1 << 3);
94    pub const SEQUENCE_CACHE: Self = Self(1 << 4);
95    pub const POLYMER_CACHE: Self = Self(1 << 5);
96    pub const ALTLOC_CACHE: Self = Self(1 << 6);
97    pub const ASSEMBLY_CACHE: Self = Self(1 << 7);
98    pub const BOND_CACHE: Self = Self(1 << 8);
99    pub const BACKBONE_GEOMETRY: Self = Self(1 << 9);
100    pub const SIDECHAIN_GEOMETRY: Self = Self(1 << 10);
101    pub const NUCLEIC_GEOMETRY: Self = Self(1 << 11);
102    pub const SECONDARY_STRUCTURE: Self = Self(1 << 12);
103    pub const CONTACT_MAP: Self = Self(1 << 13);
104    pub const GRAPH_CACHE: Self = Self(1 << 14);
105
106    #[must_use]
107    pub const fn union(self, other: Self) -> Self {
108        Self(self.0 | other.0)
109    }
110
111    #[must_use]
112    pub const fn contains(self, other: Self) -> bool {
113        (self.0 & other.0) == other.0
114    }
115}
116
117impl std::ops::BitOr for BioDerivedState {
118    type Output = Self;
119
120    fn bitor(self, rhs: Self) -> Self::Output {
121        Self(self.0 | rhs.0)
122    }
123}
124
125// ---------------------------------------------------------------------------
126// Operation classification enums
127// ---------------------------------------------------------------------------
128
129#[derive(Debug, Clone, Copy, PartialEq, Eq)]
130pub enum BioOpKind {
131    /// Does not change row identity (coordinate transforms, annotations, cache).
132    Weak,
133    /// Changes row topology or hierarchy identity (remove residues, assembly
134    /// expansion, altloc resolution, merge). Requires a mapping.
135    Strong,
136}
137
138#[derive(Debug, Clone, Copy, PartialEq, Eq)]
139pub enum BioEditKind {
140    None,
141    Local,
142    Compacting,
143    Expanding,
144    Renumbering,
145    Splitting,
146    Merging,
147    Transforming,
148}
149
150#[derive(Debug, Clone, Copy, PartialEq, Eq)]
151pub enum BioOpDomain {
152    Selection,
153    Hierarchy,
154    Coordinate,
155    Assembly,
156    Annotation,
157    Bonding,
158    Polymer,
159    ChemistryBridge,
160}
161
162#[derive(Debug, Clone, Copy, PartialEq, Eq)]
163pub enum BioParityPolicy {
164    NotApplicable,
165    GemmiWhenApplicable,
166    BiopythonWhenApplicable,
167    PdbSpecRequired,
168    RequiredNow,
169}
170
171#[derive(Debug, Clone, Copy, PartialEq, Eq)]
172pub enum MappingRequirement {
173    None,
174    Identity,
175    Required,
176}
177
178// ---------------------------------------------------------------------------
179// BioStructureOpSpec
180// ---------------------------------------------------------------------------
181
182#[derive(Debug, Clone, Copy, PartialEq, Eq)]
183pub struct BioStructureOpSpec {
184    pub method: &'static str,
185    pub impl_fn: &'static str,
186    pub domain: BioOpDomain,
187    pub kind: BioOpKind,
188    pub edit_kind: BioEditKind,
189    pub may_mutate: BioBlockSet,
190    pub auto_remap: BioBlockSet,
191    pub must_handle: BioStateSet,
192    pub needs_update: BioDerivedState,
193    pub requires_mapping: MappingRequirement,
194    pub allows_noop: bool,
195    pub support: SupportStatus,
196    pub parity: BioParityPolicy,
197    pub io_roundtrip: bool,
198}
199
200impl std::fmt::Display for BioStructureOpSpec {
201    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
202        f.write_str(self.method)
203    }
204}
205
206// ---------------------------------------------------------------------------
207// Operation outcome and errors
208// ---------------------------------------------------------------------------
209
210#[derive(Debug, Clone, PartialEq, Eq)]
211pub enum BioOpOutcome {
212    Changed,
213    NoOp { reason: &'static str },
214}
215
216#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
217pub enum BioOperationError {
218    #[error("{operation}: unsupported operation: {reason}")]
219    Unsupported {
220        operation: &'static BioStructureOpSpec,
221        reason: &'static str,
222    },
223    #[error("{operation}: invalid input: {message}")]
224    InvalidInput {
225        operation: &'static BioStructureOpSpec,
226        message: &'static str,
227    },
228    #[error("{operation}: invariant violation: {message}")]
229    InvariantViolation {
230        operation: &'static BioStructureOpSpec,
231        message: &'static str,
232    },
233}
234
235// ---------------------------------------------------------------------------
236// Mapping system
237// ---------------------------------------------------------------------------
238
239#[derive(Debug, Clone, PartialEq, Eq)]
240pub struct BioRowMapping<T: Copy> {
241    pub old_to_new: Vec<Option<T>>,
242    pub new_to_old: Vec<T>,
243}
244
245impl<T: Copy> BioRowMapping<T> {
246    #[must_use]
247    pub fn identity(len: usize, make: impl Fn(u32) -> T) -> Self {
248        let ids: Vec<T> = (0..len as u32).map(&make).collect();
249        Self {
250            old_to_new: ids.iter().copied().map(Some).collect(),
251            new_to_old: ids,
252        }
253    }
254}
255
256#[derive(Debug, Clone, PartialEq, Eq)]
257pub struct BioStructureMapping {
258    pub atoms: BioRowMapping<AtomId>,
259    pub residues: BioRowMapping<ResidueId>,
260    pub chains: BioRowMapping<ChainId>,
261    pub entities: BioRowMapping<EntityId>,
262    pub models: BioRowMapping<ModelId>,
263    pub bonds: BioRowMapping<BondId>,
264    pub assemblies: BioRowMapping<AssemblyId>,
265}
266
267// ---------------------------------------------------------------------------
268// BioOpParts — the only mutable capability object for operation bodies
269// ---------------------------------------------------------------------------
270
271#[derive(Debug, Clone, Copy, PartialEq, Eq)]
272pub struct BioOperationTrace {
273    touched_blocks: BioBlockSet,
274    remapped_blocks: BioBlockSet,
275    handled: BioStateSet,
276    cleared_cache: BioDerivedState,
277    updated_cache: BioDerivedState,
278    outcome_recorded: bool,
279}
280
281pub struct BioOpParts<'a> {
282    spec: &'static BioStructureOpSpec,
283    working: BioStructure,
284    mapping: Option<BioStructureMapping>,
285    _source: PhantomData<&'a BioStructure>,
286
287    #[cfg(feature = "op-contracts")]
288    trace: BioOperationTrace,
289}
290
291impl<'a> BioOpParts<'a> {
292    pub(crate) fn new(source: &'a BioStructure, spec: &'static BioStructureOpSpec) -> Self {
293        Self {
294            spec,
295            working: source.clone(),
296            mapping: None,
297            _source: PhantomData,
298            #[cfg(feature = "op-contracts")]
299            trace: BioOperationTrace {
300                touched_blocks: BioBlockSet::NONE,
301                remapped_blocks: BioBlockSet::NONE,
302                handled: BioStateSet::NONE,
303                cleared_cache: BioDerivedState::NONE,
304                updated_cache: BioDerivedState::NONE,
305                outcome_recorded: false,
306            },
307        }
308    }
309
310    #[must_use]
311    pub(crate) fn structure(&self) -> &BioStructure {
312        &self.working
313    }
314
315    pub(crate) fn clear_cache(&mut self, states: BioDerivedState) {
316        #[cfg(feature = "op-contracts")]
317        {
318            self.trace.cleared_cache = self.trace.cleared_cache | states;
319        }
320        let _ = states;
321    }
322
323    fn mark_handled(&mut self, states: BioStateSet) {
324        #[cfg(feature = "op-contracts")]
325        {
326            self.trace.handled = self.trace.handled.union(states);
327        }
328        let _ = states;
329    }
330
331    fn record_remapped(&mut self, blocks: BioBlockSet) {
332        #[cfg(feature = "op-contracts")]
333        {
334            self.trace.remapped_blocks = self.trace.remapped_blocks.union(blocks);
335        }
336        let _ = blocks;
337    }
338
339    pub(crate) fn record_identity_mapping(&mut self) {
340        self.mapping = Some(BioStructureMapping {
341            atoms: BioRowMapping::identity(self.working.atoms.len(), AtomId::new),
342            residues: BioRowMapping::identity(self.working.residues.len(), ResidueId::new),
343            chains: BioRowMapping::identity(self.working.chains.len(), ChainId::new),
344            entities: BioRowMapping {
345                old_to_new: vec![],
346                new_to_old: vec![],
347            },
348            models: BioRowMapping::identity(self.working.models.len(), ModelId::new),
349            bonds: BioRowMapping {
350                old_to_new: vec![],
351                new_to_old: vec![],
352            },
353            assemblies: BioRowMapping {
354                old_to_new: vec![],
355                new_to_old: vec![],
356            },
357        });
358    }
359
360    pub(crate) fn mark_hierarchy_contract_handled(&mut self) {
361        self.mark_handled(
362            BioStateSet::HIERARCHY
363                .union(BioStateSet::RESIDUE_SPANS)
364                .union(BioStateSet::CHAIN_SPANS)
365                .union(BioStateSet::MODEL_SPANS)
366                .union(BioStateSet::COORDINATE_ALIGNMENT),
367        );
368        self.record_remapped(BioBlockSet::COORDINATES);
369    }
370
371    pub(crate) fn remove_residues(
372        &mut self,
373        residues_to_remove: &[ResidueId],
374    ) -> Result<&BioStructureMapping, BioOperationError> {
375        self.assert_compacting_hierarchy_edit_allowed()?;
376
377        let mut remove_residue = vec![false; self.working.residues.len()];
378        for residue in residues_to_remove {
379            if let Some(slot) = remove_residue.get_mut(residue.index() as usize) {
380                *slot = true;
381            }
382        }
383
384        let keep_residue: Vec<bool> = remove_residue.iter().map(|remove| !remove).collect();
385        let keep_atom: Vec<bool> = self
386            .working
387            .atoms
388            .iter()
389            .map(|atom| keep_residue[atom.residue_id.index() as usize])
390            .collect();
391
392        let mut atom_old_to_new = vec![None; keep_atom.len()];
393        let mut atom_new_to_old = Vec::new();
394        for (old, keep) in keep_atom.iter().copied().enumerate() {
395            if keep {
396                let new_id = AtomId::new(atom_new_to_old.len() as u32);
397                atom_old_to_new[old] = Some(new_id);
398                atom_new_to_old.push(AtomId::new(old as u32));
399            }
400        }
401
402        let mut residue_old_to_new = vec![None; keep_residue.len()];
403        let mut residue_new_to_old = Vec::new();
404        for (old, keep) in keep_residue.iter().copied().enumerate() {
405            if keep {
406                let new_id = ResidueId::new(residue_new_to_old.len() as u32);
407                residue_old_to_new[old] = Some(new_id);
408                residue_new_to_old.push(ResidueId::new(old as u32));
409            }
410        }
411
412        let new_atoms: Vec<_> = atom_new_to_old
413            .iter()
414            .map(|old_id| {
415                let mut row = self.working.atoms[old_id.index() as usize].clone();
416                row.residue_id = residue_old_to_new[row.residue_id.index() as usize]
417                    .expect("kept atom must belong to a kept residue");
418                row
419            })
420            .collect();
421
422        let new_residues: Vec<_> = residue_new_to_old
423            .iter()
424            .map(|old_id| {
425                let old_row = &self.working.residues[old_id.index() as usize];
426                let new_start = (old_row.atom_span.start..old_row.atom_span.end())
427                    .find_map(|idx| atom_old_to_new[idx as usize].map(AtomId::index))
428                    .unwrap_or(new_atoms.len() as u32);
429                let new_len = (old_row.atom_span.start..old_row.atom_span.end())
430                    .filter(|idx| atom_old_to_new[*idx as usize].is_some())
431                    .count() as u32;
432                let mut row = old_row.clone();
433                row.atom_span = crate::bio::RowSpan::new(new_start, new_len);
434                row
435            })
436            .collect();
437
438        let new_chains: Vec<_> = self
439            .working
440            .chains
441            .iter()
442            .map(|chain| {
443                let new_start = (chain.residue_span.start..chain.residue_span.end())
444                    .find_map(|idx| residue_old_to_new[idx as usize].map(ResidueId::index))
445                    .unwrap_or(new_residues.len() as u32);
446                let new_len = (chain.residue_span.start..chain.residue_span.end())
447                    .filter(|idx| residue_old_to_new[*idx as usize].is_some())
448                    .count() as u32;
449                let mut row = chain.clone();
450                row.residue_span = crate::bio::RowSpan::new(new_start, new_len);
451                row
452            })
453            .collect();
454
455        let new_positions: Vec<_> = atom_new_to_old
456            .iter()
457            .map(|old_id| self.working.coordinates.positions[old_id.index() as usize])
458            .collect();
459
460        self.record_mutation(BioBlockSet::ATOMS);
461        self.working.atoms = new_atoms;
462        self.record_mutation(BioBlockSet::RESIDUES);
463        self.working.residues = new_residues;
464        self.record_mutation(BioBlockSet::CHAINS);
465        self.working.chains = new_chains;
466        self.record_mutation(BioBlockSet::COORDINATES);
467        self.working.coordinates.positions = new_positions;
468
469        self.mark_hierarchy_contract_handled();
470        self.clear_cache(self.spec.needs_update);
471        self.mapping = Some(BioStructureMapping {
472            atoms: BioRowMapping {
473                old_to_new: atom_old_to_new,
474                new_to_old: atom_new_to_old,
475            },
476            residues: BioRowMapping {
477                old_to_new: residue_old_to_new,
478                new_to_old: residue_new_to_old,
479            },
480            chains: BioRowMapping::identity(self.working.chains.len(), ChainId::new),
481            entities: BioRowMapping {
482                old_to_new: vec![],
483                new_to_old: vec![],
484            },
485            models: BioRowMapping::identity(self.working.models.len(), ModelId::new),
486            bonds: BioRowMapping {
487                old_to_new: vec![],
488                new_to_old: vec![],
489            },
490            assemblies: BioRowMapping {
491                old_to_new: vec![],
492                new_to_old: vec![],
493            },
494        });
495
496        Ok(self.mapping.as_ref().expect("mapping was just recorded"))
497    }
498
499    pub(crate) fn finish(
500        #[cfg_attr(not(feature = "op-contracts"), allow(unused_mut))] mut self,
501        outcome: BioOpOutcome,
502    ) -> Result<BioStructure, BioOperationError> {
503        #[cfg(feature = "op-contracts")]
504        {
505            self.trace.outcome_recorded = true;
506            let _ = outcome;
507            self.validate_contract()?;
508        }
509        #[cfg(not(feature = "op-contracts"))]
510        {
511            let _ = outcome;
512        }
513        if self.spec.requires_mapping == MappingRequirement::Required && self.mapping.is_none() {
514            return Err(BioOperationError::InvalidInput {
515                operation: self.spec,
516                message: "strong operation did not record a BioStructureMapping",
517            });
518        }
519        crate::bio_invariants::enforce_bio_structure_invariants(&self.working).map_err(
520            |message| BioOperationError::InvariantViolation {
521                operation: self.spec,
522                message,
523            },
524        )?;
525        Ok(self.working)
526    }
527
528    fn record_mutation(&mut self, block: BioBlockSet) {
529        #[cfg(feature = "op-contracts")]
530        {
531            assert!(
532                self.spec.may_mutate.contains(block),
533                "bio operation `{}` attempted to mutate a block outside its registry permissions",
534                self.spec.method
535            );
536            self.trace.touched_blocks = self.trace.touched_blocks.union(block);
537        }
538        let _ = block;
539    }
540
541    fn assert_compacting_hierarchy_edit_allowed(&self) -> Result<(), BioOperationError> {
542        if self.spec.kind != BioOpKind::Strong {
543            return Err(BioOperationError::InvalidInput {
544                operation: self.spec,
545                message: "compacting hierarchy edits require a strong operation",
546            });
547        }
548        if self.spec.edit_kind != BioEditKind::Compacting {
549            return Err(BioOperationError::InvalidInput {
550                operation: self.spec,
551                message: "operation registry does not allow compacting hierarchy edits",
552            });
553        }
554        if self.spec.requires_mapping != MappingRequirement::Required {
555            return Err(BioOperationError::InvalidInput {
556                operation: self.spec,
557                message: "compacting hierarchy edits must require a mapping",
558            });
559        }
560        Ok(())
561    }
562
563    #[cfg(feature = "op-contracts")]
564    fn validate_contract(&self) -> Result<(), BioOperationError> {
565        if !self.trace.handled.contains(self.spec.must_handle) {
566            return Err(BioOperationError::InvalidInput {
567                operation: self.spec,
568                message: "operation body did not handle every required BioStructure state",
569            });
570        }
571        let updated_or_cleared = self.trace.cleared_cache | self.trace.updated_cache;
572        if !updated_or_cleared.contains(self.spec.needs_update) {
573            return Err(BioOperationError::InvalidInput {
574                operation: self.spec,
575                message: "operation body did not clear or update every required BioStructure cache state",
576            });
577        }
578        if !self.trace.remapped_blocks.contains(self.spec.auto_remap) {
579            return Err(BioOperationError::InvalidInput {
580                operation: self.spec,
581                message: "operation did not remap every registry-required BioStructure block",
582            });
583        }
584        Ok(())
585    }
586}
587
588// ---------------------------------------------------------------------------
589// Registry tables (populated by bio_structure_ops! macro)
590// ---------------------------------------------------------------------------
591
592#[derive(Debug, Clone, Copy, PartialEq, Eq)]
593pub struct BioSupportMatrixEntry {
594    pub feature: &'static crate::FeatureSpec,
595    pub operation: &'static BioStructureOpSpec,
596}
597
598#[derive(Debug, Clone, Copy, PartialEq, Eq)]
599pub struct BioOperationInvariantEntry {
600    pub operation: &'static BioStructureOpSpec,
601    pub profile: &'static str,
602}
603
604#[derive(Debug, Clone, Copy, PartialEq, Eq)]
605pub struct BioParityMatrixEntry {
606    pub operation: &'static BioStructureOpSpec,
607    pub profile: &'static str,
608}
609
610// ---------------------------------------------------------------------------
611// Operation registry
612// ---------------------------------------------------------------------------
613
614use cosmolkit_macros::bio_op_body;
615use cosmolkit_macros::bio_structure_ops;
616
617bio_structure_ops! {
618    op remove_waters() {
619        method: without_waters,
620        impl_fn: remove_waters_impl,
621        domain: selection,
622        kind: strong,
623        edit_kind: compacting,
624        may_mutate: [atoms, residues, chains, models, coordinates],
625        auto_remap: [coordinates],
626        must_handle: [hierarchy, residue_spans, chain_spans, model_spans, coordinate_alignment],
627        needs_update: [atom_index, residue_index, chain_index],
628        requires_mapping: required,
629        allows_noop: true,
630        feature: BIO_SELECTION_FEATURE,
631        parity: not_applicable,
632        io_roundtrip: false,
633        invariant_profile: "strong_bio_hierarchy",
634    }
635}
636
637#[bio_op_body(remove_waters, parts)]
638fn remove_waters_impl() -> Result<BioOpOutcome, BioOperationError> {
639    use crate::bio::ResidueKind;
640
641    let water_residue_ids: Vec<ResidueId> = parts
642        .structure()
643        .residues
644        .iter()
645        .enumerate()
646        .filter(|(_, r)| r.kind == ResidueKind::Water)
647        .map(|(index, _)| ResidueId::new(index as u32))
648        .collect();
649
650    if water_residue_ids.is_empty() {
651        parts.record_identity_mapping();
652        parts.mark_hierarchy_contract_handled();
653        parts.clear_cache(BIO_REMOVE_WATERS_SPEC.needs_update);
654        return Ok(BioOpOutcome::NoOp {
655            reason: "no water residues found",
656        });
657    }
658
659    parts.remove_residues(&water_residue_ids)?;
660    Ok(BioOpOutcome::Changed)
661}
662
663#[cfg(test)]
664mod tests {
665    use super::*;
666    use crate::bio::*;
667
668    const TEST_WEAK_COMPACT_SPEC: BioStructureOpSpec = BioStructureOpSpec {
669        method: "test_weak_compact",
670        impl_fn: "test_weak_compact_impl",
671        domain: BioOpDomain::Selection,
672        kind: BioOpKind::Weak,
673        edit_kind: BioEditKind::Compacting,
674        may_mutate: BioBlockSet::ATOMS,
675        auto_remap: BioBlockSet::NONE,
676        must_handle: BioStateSet::NONE,
677        needs_update: BioDerivedState::NONE,
678        requires_mapping: MappingRequirement::Required,
679        allows_noop: true,
680        support: SupportStatus::Experimental,
681        parity: BioParityPolicy::NotApplicable,
682        io_roundtrip: false,
683    };
684
685    #[cfg(feature = "op-contracts")]
686    const TEST_UNAUTHORIZED_REMOVE_RESIDUES_SPEC: BioStructureOpSpec = BioStructureOpSpec {
687        method: "test_unauthorized_remove_residues",
688        impl_fn: "test_unauthorized_remove_residues_impl",
689        domain: BioOpDomain::Selection,
690        kind: BioOpKind::Strong,
691        edit_kind: BioEditKind::Compacting,
692        may_mutate: BioBlockSet::NONE,
693        auto_remap: BioBlockSet::NONE,
694        must_handle: BioStateSet::NONE,
695        needs_update: BioDerivedState::NONE,
696        requires_mapping: MappingRequirement::Required,
697        allows_noop: true,
698        support: SupportStatus::Experimental,
699        parity: BioParityPolicy::NotApplicable,
700        io_roundtrip: false,
701    };
702
703    fn make_structure_with_waters() -> BioStructure {
704        // Model 0 → Chain 0 → Residue 0 (ALA, atom 0), Residue 1 (HOH, atom 1)
705        let mut s = BioStructure::new();
706        s.models.push(ModelRow {
707            chain_span: RowSpan::new(0, 1),
708            source_model_number: Some(1),
709        });
710        s.chains.push(ChainRow {
711            model_id: ModelId::new(0),
712            entity_id: None,
713            residue_span: RowSpan::new(0, 2),
714            kind: ChainKind::Mixed,
715            source: ChainSourceIds {
716                auth_chain_id: None,
717                label_asym_id: None,
718            },
719        });
720        s.residues.push(ResidueRow {
721            chain_id: ChainId::new(0),
722            atom_span: RowSpan::new(0, 1),
723            name: ResidueName([b'A', b'L', b'A', 0], 3),
724            kind: ResidueKind::AminoAcid,
725            entity_kind: EntityKind::Unknown,
726            source: ResidueSourceIds {
727                seq_id: None,
728                label_seq_id: None,
729                segment_id: None,
730                subchain_id: None,
731                label_entity_id: None,
732            },
733            het_flag: None,
734            sifts_unp: None,
735        });
736        s.residues.push(ResidueRow {
737            chain_id: ChainId::new(0),
738            atom_span: RowSpan::new(1, 1),
739            name: ResidueName([b'H', b'O', b'H', 0], 3),
740            kind: ResidueKind::Water,
741            entity_kind: EntityKind::Unknown,
742            source: ResidueSourceIds {
743                seq_id: None,
744                label_seq_id: None,
745                segment_id: None,
746                subchain_id: None,
747                label_entity_id: None,
748            },
749            het_flag: None,
750            sifts_unp: None,
751        });
752        s.atoms.push(AtomRow {
753            residue_id: ResidueId::new(0),
754            name: AtomName([b' ', b'C', b'A', b' ']),
755            element: crate::Element::C,
756            altloc: None,
757            occupancy: None,
758            b_iso: None,
759            formal_charge: None,
760            anisou: None,
761            calc_flag: BioCalcFlag::NotSet,
762            tls_group_id: None,
763            fraction: None,
764            source: AtomSourceIds { serial: None },
765        });
766        s.atoms.push(AtomRow {
767            residue_id: ResidueId::new(1),
768            name: AtomName([b' ', b'O', b' ', b' ']),
769            element: crate::Element::O,
770            altloc: None,
771            occupancy: None,
772            b_iso: None,
773            formal_charge: None,
774            anisou: None,
775            calc_flag: BioCalcFlag::NotSet,
776            tls_group_id: None,
777            fraction: None,
778            source: AtomSourceIds { serial: None },
779        });
780        s.coordinates.positions = vec![[1.0, 0.0, 0.0], [5.0, 0.0, 0.0]];
781        s
782    }
783
784    #[test]
785    fn registered_bio_ops_have_matrix_entries() {
786        assert_eq!(BIO_STRUCTURE_OPS.len(), 1);
787        for operation in BIO_STRUCTURE_OPS {
788            assert!(
789                BIO_SUPPORT_MATRIX
790                    .iter()
791                    .any(|entry| std::ptr::eq(entry.operation, *operation)),
792                "missing bio support matrix entry for {}",
793                operation.method
794            );
795            assert!(
796                BIO_OPERATION_INVARIANT_MATRIX
797                    .iter()
798                    .any(|entry| std::ptr::eq(entry.operation, *operation)),
799                "missing bio invariant matrix entry for {}",
800                operation.method
801            );
802        }
803    }
804
805    #[test]
806    fn remove_waters_removes_water_residue_and_atom() {
807        let s = make_structure_with_waters();
808        let result = s.without_waters().expect("remove_waters should succeed");
809
810        assert_eq!(result.num_atoms(), 1);
811        assert_eq!(result.num_residues(), 1);
812        assert_eq!(result.residues[0].kind, ResidueKind::AminoAcid);
813        assert_eq!(result.coordinates.positions, vec![[1.0, 0.0, 0.0]]);
814    }
815
816    #[test]
817    fn remove_waters_is_noop_on_structure_without_waters() {
818        let mut s = BioStructure::new();
819        s.models.push(ModelRow {
820            chain_span: RowSpan::new(0, 0),
821            source_model_number: None,
822        });
823        let result = s.without_waters().expect("noop should succeed");
824        assert_eq!(result.num_atoms(), 0);
825    }
826
827    #[test]
828    fn remove_waters_preserves_source_invariants() {
829        let s = make_structure_with_waters();
830        let result = s.without_waters().unwrap();
831        crate::bio_invariants::enforce_bio_structure_invariants(&result)
832            .expect("result must satisfy invariants");
833    }
834
835    #[test]
836    fn remove_residues_rejects_weak_operation_specs() {
837        let s = BioStructure::new();
838        let mut parts = BioOpParts::new(&s, &TEST_WEAK_COMPACT_SPEC);
839        let err = parts
840            .remove_residues(&[])
841            .expect_err("weak operation must not compact hierarchy");
842        assert!(matches!(err, BioOperationError::InvalidInput { .. }));
843    }
844
845    #[cfg(feature = "op-contracts")]
846    #[test]
847    #[should_panic(expected = "attempted to mutate a block outside its registry permissions")]
848    fn remove_residues_panics_when_registry_does_not_allow_mutation() {
849        let s = make_structure_with_waters();
850        let mut parts = BioOpParts::new(&s, &TEST_UNAUTHORIZED_REMOVE_RESIDUES_SPEC);
851        let water = ResidueId::new(1);
852        let _ = parts.remove_residues(&[water]);
853    }
854}