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 support: SupportStatus,
195    pub parity: BioParityPolicy,
196    pub io_roundtrip: bool,
197}
198
199impl std::fmt::Display for BioStructureOpSpec {
200    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
201        f.write_str(self.method)
202    }
203}
204
205// ---------------------------------------------------------------------------
206// Operation errors
207// ---------------------------------------------------------------------------
208
209#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
210pub enum BioOperationError {
211    #[error("{operation}: unsupported operation: {reason}")]
212    Unsupported {
213        operation: &'static BioStructureOpSpec,
214        reason: &'static str,
215    },
216    #[error("{operation}: invalid input: {message}")]
217    InvalidInput {
218        operation: &'static BioStructureOpSpec,
219        message: &'static str,
220    },
221    #[error("{operation}: invariant violation: {message}")]
222    InvariantViolation {
223        operation: &'static BioStructureOpSpec,
224        message: &'static str,
225    },
226}
227
228// ---------------------------------------------------------------------------
229// Mapping system
230// ---------------------------------------------------------------------------
231
232#[derive(Debug, Clone, PartialEq, Eq)]
233pub struct BioRowMapping<T: Copy> {
234    pub old_to_new: Vec<Option<T>>,
235    pub new_to_old: Vec<T>,
236}
237
238impl<T: Copy> BioRowMapping<T> {
239    #[must_use]
240    pub fn identity(len: usize, make: impl Fn(u32) -> T) -> Self {
241        let ids: Vec<T> = (0..len as u32).map(&make).collect();
242        Self {
243            old_to_new: ids.iter().copied().map(Some).collect(),
244            new_to_old: ids,
245        }
246    }
247}
248
249#[derive(Debug, Clone, PartialEq, Eq)]
250pub struct BioStructureMapping {
251    pub atoms: BioRowMapping<AtomId>,
252    pub residues: BioRowMapping<ResidueId>,
253    pub chains: BioRowMapping<ChainId>,
254    pub entities: BioRowMapping<EntityId>,
255    pub models: BioRowMapping<ModelId>,
256    pub bonds: BioRowMapping<BondId>,
257    pub assemblies: BioRowMapping<AssemblyId>,
258}
259
260// ---------------------------------------------------------------------------
261// BioOpParts — the only mutable capability object for operation bodies
262// ---------------------------------------------------------------------------
263
264#[derive(Debug, Clone, Copy, PartialEq, Eq)]
265pub struct BioOperationTrace {
266    touched_blocks: BioBlockSet,
267    remapped_blocks: BioBlockSet,
268    handled: BioStateSet,
269    cleared_cache: BioDerivedState,
270    updated_cache: BioDerivedState,
271}
272
273pub struct BioOpParts<'a> {
274    spec: &'static BioStructureOpSpec,
275    working: BioStructure,
276    mapping: Option<BioStructureMapping>,
277    _source: PhantomData<&'a BioStructure>,
278
279    #[cfg(feature = "op-contracts")]
280    trace: BioOperationTrace,
281}
282
283impl<'a> BioOpParts<'a> {
284    pub(crate) fn new(source: &'a BioStructure, spec: &'static BioStructureOpSpec) -> Self {
285        Self {
286            spec,
287            working: source.clone(),
288            mapping: None,
289            _source: PhantomData,
290            #[cfg(feature = "op-contracts")]
291            trace: BioOperationTrace {
292                touched_blocks: BioBlockSet::NONE,
293                remapped_blocks: BioBlockSet::NONE,
294                handled: BioStateSet::NONE,
295                cleared_cache: BioDerivedState::NONE,
296                updated_cache: BioDerivedState::NONE,
297            },
298        }
299    }
300
301    #[must_use]
302    pub(crate) fn structure(&self) -> &BioStructure {
303        &self.working
304    }
305
306    pub(crate) fn clear_cache(&mut self, states: BioDerivedState) {
307        #[cfg(feature = "op-contracts")]
308        {
309            self.trace.cleared_cache = self.trace.cleared_cache | states;
310        }
311        let _ = states;
312    }
313
314    fn mark_handled(&mut self, states: BioStateSet) {
315        #[cfg(feature = "op-contracts")]
316        {
317            self.trace.handled = self.trace.handled.union(states);
318        }
319        let _ = states;
320    }
321
322    fn record_remapped(&mut self, blocks: BioBlockSet) {
323        #[cfg(feature = "op-contracts")]
324        {
325            self.trace.remapped_blocks = self.trace.remapped_blocks.union(blocks);
326        }
327        let _ = blocks;
328    }
329
330    pub(crate) fn record_identity_mapping(&mut self) {
331        self.mapping = Some(BioStructureMapping {
332            atoms: BioRowMapping::identity(self.working.atoms.len(), AtomId::new),
333            residues: BioRowMapping::identity(self.working.residues.len(), ResidueId::new),
334            chains: BioRowMapping::identity(self.working.chains.len(), ChainId::new),
335            entities: BioRowMapping {
336                old_to_new: vec![],
337                new_to_old: vec![],
338            },
339            models: BioRowMapping::identity(self.working.models.len(), ModelId::new),
340            bonds: BioRowMapping {
341                old_to_new: vec![],
342                new_to_old: vec![],
343            },
344            assemblies: BioRowMapping {
345                old_to_new: vec![],
346                new_to_old: vec![],
347            },
348        });
349    }
350
351    pub(crate) fn mark_hierarchy_contract_handled(&mut self) {
352        self.mark_handled(
353            BioStateSet::HIERARCHY
354                .union(BioStateSet::RESIDUE_SPANS)
355                .union(BioStateSet::CHAIN_SPANS)
356                .union(BioStateSet::MODEL_SPANS)
357                .union(BioStateSet::COORDINATE_ALIGNMENT),
358        );
359        self.record_remapped(BioBlockSet::COORDINATES);
360    }
361
362    pub(crate) fn remove_residues(
363        &mut self,
364        residues_to_remove: &[ResidueId],
365    ) -> Result<&BioStructureMapping, BioOperationError> {
366        self.assert_compacting_hierarchy_edit_allowed()?;
367
368        let mut remove_residue = vec![false; self.working.residues.len()];
369        for residue in residues_to_remove {
370            if let Some(slot) = remove_residue.get_mut(residue.index() as usize) {
371                *slot = true;
372            }
373        }
374
375        let keep_residue: Vec<bool> = remove_residue.iter().map(|remove| !remove).collect();
376        let keep_atom: Vec<bool> = self
377            .working
378            .atoms
379            .iter()
380            .map(|atom| keep_residue[atom.residue_id.index() as usize])
381            .collect();
382
383        let mut atom_old_to_new = vec![None; keep_atom.len()];
384        let mut atom_new_to_old = Vec::new();
385        for (old, keep) in keep_atom.iter().copied().enumerate() {
386            if keep {
387                let new_id = AtomId::new(atom_new_to_old.len() as u32);
388                atom_old_to_new[old] = Some(new_id);
389                atom_new_to_old.push(AtomId::new(old as u32));
390            }
391        }
392
393        let mut residue_old_to_new = vec![None; keep_residue.len()];
394        let mut residue_new_to_old = Vec::new();
395        for (old, keep) in keep_residue.iter().copied().enumerate() {
396            if keep {
397                let new_id = ResidueId::new(residue_new_to_old.len() as u32);
398                residue_old_to_new[old] = Some(new_id);
399                residue_new_to_old.push(ResidueId::new(old as u32));
400            }
401        }
402
403        let new_atoms: Vec<_> = atom_new_to_old
404            .iter()
405            .map(|old_id| {
406                let mut row = self.working.atoms[old_id.index() as usize].clone();
407                row.residue_id = residue_old_to_new[row.residue_id.index() as usize]
408                    .expect("kept atom must belong to a kept residue");
409                row
410            })
411            .collect();
412
413        let new_residues: Vec<_> = residue_new_to_old
414            .iter()
415            .map(|old_id| {
416                let old_row = &self.working.residues[old_id.index() as usize];
417                let new_start = (old_row.atom_span.start..old_row.atom_span.end())
418                    .find_map(|idx| atom_old_to_new[idx as usize].map(AtomId::index))
419                    .unwrap_or(new_atoms.len() as u32);
420                let new_len = (old_row.atom_span.start..old_row.atom_span.end())
421                    .filter(|idx| atom_old_to_new[*idx as usize].is_some())
422                    .count() as u32;
423                let mut row = old_row.clone();
424                row.atom_span = crate::bio::RowSpan::new(new_start, new_len);
425                row
426            })
427            .collect();
428
429        let new_chains: Vec<_> = self
430            .working
431            .chains
432            .iter()
433            .map(|chain| {
434                let new_start = (chain.residue_span.start..chain.residue_span.end())
435                    .find_map(|idx| residue_old_to_new[idx as usize].map(ResidueId::index))
436                    .unwrap_or(new_residues.len() as u32);
437                let new_len = (chain.residue_span.start..chain.residue_span.end())
438                    .filter(|idx| residue_old_to_new[*idx as usize].is_some())
439                    .count() as u32;
440                let mut row = chain.clone();
441                row.residue_span = crate::bio::RowSpan::new(new_start, new_len);
442                row
443            })
444            .collect();
445
446        let new_positions: Vec<_> = atom_new_to_old
447            .iter()
448            .map(|old_id| self.working.coordinates.positions[old_id.index() as usize])
449            .collect();
450
451        self.record_mutation(BioBlockSet::ATOMS);
452        self.working.atoms = new_atoms;
453        self.record_mutation(BioBlockSet::RESIDUES);
454        self.working.residues = new_residues;
455        self.record_mutation(BioBlockSet::CHAINS);
456        self.working.chains = new_chains;
457        self.record_mutation(BioBlockSet::COORDINATES);
458        self.working.coordinates.positions = new_positions;
459
460        self.mark_hierarchy_contract_handled();
461        self.clear_cache(self.spec.needs_update);
462        self.mapping = Some(BioStructureMapping {
463            atoms: BioRowMapping {
464                old_to_new: atom_old_to_new,
465                new_to_old: atom_new_to_old,
466            },
467            residues: BioRowMapping {
468                old_to_new: residue_old_to_new,
469                new_to_old: residue_new_to_old,
470            },
471            chains: BioRowMapping::identity(self.working.chains.len(), ChainId::new),
472            entities: BioRowMapping {
473                old_to_new: vec![],
474                new_to_old: vec![],
475            },
476            models: BioRowMapping::identity(self.working.models.len(), ModelId::new),
477            bonds: BioRowMapping {
478                old_to_new: vec![],
479                new_to_old: vec![],
480            },
481            assemblies: BioRowMapping {
482                old_to_new: vec![],
483                new_to_old: vec![],
484            },
485        });
486
487        Ok(self.mapping.as_ref().expect("mapping was just recorded"))
488    }
489
490    pub(crate) fn finish(self) -> Result<BioStructure, BioOperationError> {
491        #[cfg(feature = "op-contracts")]
492        {
493            self.validate_contract()?;
494        }
495        if self.spec.requires_mapping == MappingRequirement::Required && self.mapping.is_none() {
496            return Err(BioOperationError::InvalidInput {
497                operation: self.spec,
498                message: "strong operation did not record a BioStructureMapping",
499            });
500        }
501        crate::bio_invariants::enforce_bio_structure_invariants(&self.working).map_err(
502            |message| BioOperationError::InvariantViolation {
503                operation: self.spec,
504                message,
505            },
506        )?;
507        Ok(self.working)
508    }
509
510    fn record_mutation(&mut self, block: BioBlockSet) {
511        #[cfg(feature = "op-contracts")]
512        {
513            assert!(
514                self.spec.may_mutate.contains(block),
515                "bio operation `{}` attempted to mutate a block outside its registry permissions",
516                self.spec.method
517            );
518            self.trace.touched_blocks = self.trace.touched_blocks.union(block);
519        }
520        let _ = block;
521    }
522
523    fn assert_compacting_hierarchy_edit_allowed(&self) -> Result<(), BioOperationError> {
524        if self.spec.kind != BioOpKind::Strong {
525            return Err(BioOperationError::InvalidInput {
526                operation: self.spec,
527                message: "compacting hierarchy edits require a strong operation",
528            });
529        }
530        if self.spec.edit_kind != BioEditKind::Compacting {
531            return Err(BioOperationError::InvalidInput {
532                operation: self.spec,
533                message: "operation registry does not allow compacting hierarchy edits",
534            });
535        }
536        if self.spec.requires_mapping != MappingRequirement::Required {
537            return Err(BioOperationError::InvalidInput {
538                operation: self.spec,
539                message: "compacting hierarchy edits must require a mapping",
540            });
541        }
542        Ok(())
543    }
544
545    #[cfg(feature = "op-contracts")]
546    fn validate_contract(&self) -> Result<(), BioOperationError> {
547        if !self.trace.handled.contains(self.spec.must_handle) {
548            return Err(BioOperationError::InvalidInput {
549                operation: self.spec,
550                message: "operation body did not handle every required BioStructure state",
551            });
552        }
553        let updated_or_cleared = self.trace.cleared_cache | self.trace.updated_cache;
554        if !updated_or_cleared.contains(self.spec.needs_update) {
555            return Err(BioOperationError::InvalidInput {
556                operation: self.spec,
557                message: "operation body did not clear or update every required BioStructure cache state",
558            });
559        }
560        if !self.trace.remapped_blocks.contains(self.spec.auto_remap) {
561            return Err(BioOperationError::InvalidInput {
562                operation: self.spec,
563                message: "operation did not remap every registry-required BioStructure block",
564            });
565        }
566        Ok(())
567    }
568}
569
570// ---------------------------------------------------------------------------
571// Registry tables (populated by bio_structure_ops! macro)
572// ---------------------------------------------------------------------------
573
574#[derive(Debug, Clone, Copy, PartialEq, Eq)]
575pub struct BioSupportMatrixEntry {
576    pub feature: &'static crate::FeatureSpec,
577    pub operation: &'static BioStructureOpSpec,
578}
579
580#[derive(Debug, Clone, Copy, PartialEq, Eq)]
581pub struct BioOperationInvariantEntry {
582    pub operation: &'static BioStructureOpSpec,
583    pub profile: &'static str,
584}
585
586#[derive(Debug, Clone, Copy, PartialEq, Eq)]
587pub struct BioParityMatrixEntry {
588    pub operation: &'static BioStructureOpSpec,
589    pub profile: &'static str,
590}
591
592// ---------------------------------------------------------------------------
593// Operation registry
594// ---------------------------------------------------------------------------
595
596use cosmolkit_macros::bio_op_body;
597use cosmolkit_macros::bio_structure_ops;
598
599bio_structure_ops! {
600    op remove_waters() {
601        method: without_waters,
602        impl_fn: remove_waters_impl,
603        domain: selection,
604        kind: strong,
605        edit_kind: compacting,
606        may_mutate: [atoms, residues, chains, models, coordinates],
607        auto_remap: [coordinates],
608        must_handle: [hierarchy, residue_spans, chain_spans, model_spans, coordinate_alignment],
609        needs_update: [atom_index, residue_index, chain_index],
610        requires_mapping: required,
611        feature: BIO_SELECTION_FEATURE,
612        parity: not_applicable,
613        io_roundtrip: false,
614        invariant_profile: "strong_bio_hierarchy",
615    }
616}
617
618#[bio_op_body(remove_waters, parts)]
619fn remove_waters_impl() -> Result<(), BioOperationError> {
620    use crate::bio::ResidueKind;
621
622    let water_residue_ids: Vec<ResidueId> = parts
623        .structure()
624        .residues
625        .iter()
626        .enumerate()
627        .filter(|(_, r)| r.kind == ResidueKind::Water)
628        .map(|(index, _)| ResidueId::new(index as u32))
629        .collect();
630
631    if water_residue_ids.is_empty() {
632        parts.record_identity_mapping();
633        parts.mark_hierarchy_contract_handled();
634        parts.clear_cache(BIO_REMOVE_WATERS_SPEC.needs_update);
635        return Ok(());
636    }
637
638    parts.remove_residues(&water_residue_ids)?;
639    Ok(())
640}
641
642#[cfg(test)]
643mod tests {
644    use super::*;
645    use crate::bio::*;
646
647    const TEST_WEAK_COMPACT_SPEC: BioStructureOpSpec = BioStructureOpSpec {
648        method: "test_weak_compact",
649        impl_fn: "test_weak_compact_impl",
650        domain: BioOpDomain::Selection,
651        kind: BioOpKind::Weak,
652        edit_kind: BioEditKind::Compacting,
653        may_mutate: BioBlockSet::ATOMS,
654        auto_remap: BioBlockSet::NONE,
655        must_handle: BioStateSet::NONE,
656        needs_update: BioDerivedState::NONE,
657        requires_mapping: MappingRequirement::Required,
658        support: SupportStatus::Experimental,
659        parity: BioParityPolicy::NotApplicable,
660        io_roundtrip: false,
661    };
662
663    #[cfg(feature = "op-contracts")]
664    const TEST_UNAUTHORIZED_REMOVE_RESIDUES_SPEC: BioStructureOpSpec = BioStructureOpSpec {
665        method: "test_unauthorized_remove_residues",
666        impl_fn: "test_unauthorized_remove_residues_impl",
667        domain: BioOpDomain::Selection,
668        kind: BioOpKind::Strong,
669        edit_kind: BioEditKind::Compacting,
670        may_mutate: BioBlockSet::NONE,
671        auto_remap: BioBlockSet::NONE,
672        must_handle: BioStateSet::NONE,
673        needs_update: BioDerivedState::NONE,
674        requires_mapping: MappingRequirement::Required,
675        support: SupportStatus::Experimental,
676        parity: BioParityPolicy::NotApplicable,
677        io_roundtrip: false,
678    };
679
680    fn make_structure_with_waters() -> BioStructure {
681        // Model 0 → Chain 0 → Residue 0 (ALA, atom 0), Residue 1 (HOH, atom 1)
682        let mut s = BioStructure::new();
683        s.models.push(ModelRow {
684            chain_span: RowSpan::new(0, 1),
685            source_model_number: Some(1),
686        });
687        s.chains.push(ChainRow {
688            model_id: ModelId::new(0),
689            entity_id: None,
690            residue_span: RowSpan::new(0, 2),
691            kind: ChainKind::Mixed,
692            source: ChainSourceIds {
693                auth_chain_id: None,
694                label_asym_id: None,
695            },
696        });
697        s.residues.push(ResidueRow {
698            chain_id: ChainId::new(0),
699            atom_span: RowSpan::new(0, 1),
700            name: ResidueName([b'A', b'L', b'A', 0], 3),
701            kind: ResidueKind::AminoAcid,
702            entity_kind: EntityKind::Unknown,
703            source: ResidueSourceIds {
704                seq_id: None,
705                label_seq_id: None,
706                segment_id: None,
707                subchain_id: None,
708                label_entity_id: None,
709            },
710            het_flag: None,
711            sifts_unp: None,
712        });
713        s.residues.push(ResidueRow {
714            chain_id: ChainId::new(0),
715            atom_span: RowSpan::new(1, 1),
716            name: ResidueName([b'H', b'O', b'H', 0], 3),
717            kind: ResidueKind::Water,
718            entity_kind: EntityKind::Unknown,
719            source: ResidueSourceIds {
720                seq_id: None,
721                label_seq_id: None,
722                segment_id: None,
723                subchain_id: None,
724                label_entity_id: None,
725            },
726            het_flag: None,
727            sifts_unp: None,
728        });
729        s.atoms.push(AtomRow {
730            residue_id: ResidueId::new(0),
731            name: AtomName([b' ', b'C', b'A', b' ']),
732            element: crate::Element::C,
733            altloc: None,
734            occupancy: None,
735            b_iso: None,
736            formal_charge: None,
737            anisou: None,
738            calc_flag: BioCalcFlag::NotSet,
739            tls_group_id: None,
740            fraction: None,
741            source: AtomSourceIds { serial: None },
742        });
743        s.atoms.push(AtomRow {
744            residue_id: ResidueId::new(1),
745            name: AtomName([b' ', b'O', b' ', b' ']),
746            element: crate::Element::O,
747            altloc: None,
748            occupancy: None,
749            b_iso: None,
750            formal_charge: None,
751            anisou: None,
752            calc_flag: BioCalcFlag::NotSet,
753            tls_group_id: None,
754            fraction: None,
755            source: AtomSourceIds { serial: None },
756        });
757        s.coordinates.positions = vec![[1.0, 0.0, 0.0], [5.0, 0.0, 0.0]];
758        s
759    }
760
761    #[test]
762    fn registered_bio_ops_have_matrix_entries() {
763        assert_eq!(BIO_STRUCTURE_OPS.len(), 1);
764        for operation in BIO_STRUCTURE_OPS {
765            assert!(
766                BIO_SUPPORT_MATRIX
767                    .iter()
768                    .any(|entry| std::ptr::eq(entry.operation, *operation)),
769                "missing bio support matrix entry for {}",
770                operation.method
771            );
772            assert!(
773                BIO_OPERATION_INVARIANT_MATRIX
774                    .iter()
775                    .any(|entry| std::ptr::eq(entry.operation, *operation)),
776                "missing bio invariant matrix entry for {}",
777                operation.method
778            );
779        }
780    }
781
782    #[test]
783    fn remove_waters_removes_water_residue_and_atom() {
784        let s = make_structure_with_waters();
785        let result = s.without_waters().expect("remove_waters should succeed");
786
787        assert_eq!(result.num_atoms(), 1);
788        assert_eq!(result.num_residues(), 1);
789        assert_eq!(result.residues[0].kind, ResidueKind::AminoAcid);
790        assert_eq!(result.coordinates.positions, vec![[1.0, 0.0, 0.0]]);
791    }
792
793    #[test]
794    fn remove_waters_is_noop_on_structure_without_waters() {
795        let mut s = BioStructure::new();
796        s.models.push(ModelRow {
797            chain_span: RowSpan::new(0, 0),
798            source_model_number: None,
799        });
800        let result = s.without_waters().expect("noop should succeed");
801        assert_eq!(result.num_atoms(), 0);
802    }
803
804    #[test]
805    fn remove_waters_preserves_source_invariants() {
806        let s = make_structure_with_waters();
807        let result = s.without_waters().unwrap();
808        crate::bio_invariants::enforce_bio_structure_invariants(&result)
809            .expect("result must satisfy invariants");
810    }
811
812    #[test]
813    fn remove_residues_rejects_weak_operation_specs() {
814        let s = BioStructure::new();
815        let mut parts = BioOpParts::new(&s, &TEST_WEAK_COMPACT_SPEC);
816        let err = parts
817            .remove_residues(&[])
818            .expect_err("weak operation must not compact hierarchy");
819        assert!(matches!(err, BioOperationError::InvalidInput { .. }));
820    }
821
822    #[cfg(feature = "op-contracts")]
823    #[test]
824    #[should_panic(expected = "attempted to mutate a block outside its registry permissions")]
825    fn remove_residues_panics_when_registry_does_not_allow_mutation() {
826        let s = make_structure_with_waters();
827        let mut parts = BioOpParts::new(&s, &TEST_UNAUTHORIZED_REMOVE_RESIDUES_SPEC);
828        let water = ResidueId::new(1);
829        let _ = parts.remove_residues(&[water]);
830    }
831}