use crate::FingerprintError;
use super::reaccs::MoleculeState;
use super::rings::{combine_rings, ring_list};
const MAX_NEIGHBOURS: usize = 20;
const NONE: i32 = 0;
const DOUBLET: i32 = 2;
const ZERO_COUNT: i32 = 1;
const SINGLE: i32 = 1;
const DOUBLE: i32 = 2;
const TRIPLE: i32 = 3;
const AROMATIC: i32 = 4;
const SUB_AS_IS: i32 = -2;
#[derive(Debug, Clone, Copy)]
struct ValenceEntry {
atom_type: &'static str,
from_valence: i32,
to_valence: i32,
step_valence: i32,
lone_pairs: i32,
pair_deficit: i32,
}
const VALENCE_TABLE: &[ValenceEntry] = &[
valence("C", 4, 4, 1, 0, 0),
valence("H", -1, 1, 2, 0, 0),
valence("N", 3, 5, 2, 1, 0),
valence("O", 2, 2, 1, 2, 0),
valence("Cl", 1, 7, 2, 0, 0),
valence("P", 3, 5, 2, 1, 0),
valence("S", 2, 6, 2, 2, 0),
valence("F", 1, 1, 1, 1, 0),
valence("H", -1, 1, 2, 0, 0),
valence("Li", -1, 1, 2, 0, 0),
valence("Na", -1, 1, 2, 0, 0),
valence("K", -1, 1, 2, 0, 0),
valence("Rb", -1, 1, 2, 0, 0),
valence("Cs", -1, 1, 2, 0, 0),
valence("Be", -2, 2, 2, 0, 0),
valence("Mg", -2, 2, 2, 0, 0),
valence("B", 3, 3, 1, 0, 1),
valence("Al", -3, 3, 2, 0, 1),
valence("Ga", -3, 3, 2, 0, 1),
valence("In", -3, 3, 2, 0, 1),
valence("Tl", -3, 3, 2, 0, 0),
valence("Si", 4, 4, 1, 0, 0),
valence("As", 3, 5, 2, 0, 0),
valence("Sb", 3, 5, 2, 0, 0),
valence("Bi", 3, 5, 2, 0, 0),
valence("Se", 2, 6, 2, 0, 0),
valence("Te", 2, 6, 2, 0, 0),
valence("La", -3, 3, 2, 0, 1),
valence("Ce", -3, 3, 2, 0, 1),
valence("Pr", -3, 3, 2, 0, 1),
valence("Nd", -3, 3, 2, 0, 1),
valence("Pm", -3, 3, 2, 0, 1),
valence("Sm", -3, 3, 2, 0, 1),
valence("Eu", -3, 3, 2, 0, 1),
valence("Gd", -3, 3, 2, 0, 1),
valence("Tb", -3, 3, 2, 0, 1),
valence("Dy", -3, 3, 2, 0, 1),
valence("Ho", -3, 3, 2, 0, 1),
valence("Er", -3, 3, 2, 0, 1),
valence("Tm", -3, 3, 2, 0, 1),
valence("Yb", -3, 3, 2, 0, 1),
valence("Lu", -3, 3, 2, 0, 1),
valence("Br", 1, 7, 2, 0, 0),
valence("I", 1, 7, 2, 0, 0),
];
const fn valence(
atom_type: &'static str,
from_valence: i32,
to_valence: i32,
step_valence: i32,
lone_pairs: i32,
pair_deficit: i32,
) -> ValenceEntry {
ValenceEntry {
atom_type,
from_valence,
to_valence,
step_valence,
lone_pairs,
pair_deficit,
}
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(super) struct Neighbourhood {
n_ligands: usize,
atoms: [usize; MAX_NEIGHBOURS],
bonds: [usize; MAX_NEIGHBOURS],
}
impl Default for Neighbourhood {
fn default() -> Self {
Self {
n_ligands: 0,
atoms: [0; MAX_NEIGHBOURS],
bonds: [0; MAX_NEIGHBOURS],
}
}
}
impl Neighbourhood {
pub(super) fn atoms(&self) -> &[usize] {
&self.atoms[..self.n_ligands]
}
pub(super) fn bonds(&self) -> &[usize] {
&self.bonds[..self.n_ligands]
}
fn push(&mut self, atom: usize, bond: usize, owner: usize) -> Result<(), FingerprintError> {
if self.n_ligands == MAX_NEIGHBOURS {
return Err(FingerprintError::AvalonConversion {
reason: format!(
"Avalon neighbourhood capacity exceeded at atom {}",
owner + 1
),
});
}
self.atoms[self.n_ligands] = atom;
self.bonds[self.n_ligands] = bond;
self.n_ligands += 1;
Ok(())
}
}
pub(super) fn setup_neighbourhood(
molecule: &MoleculeState,
nlimit: usize,
) -> Result<Vec<Neighbourhood>, FingerprintError> {
let mut neighbours = vec![Neighbourhood::default(); molecule.atoms.len()];
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
let at0 = source_atom_index(bond.atoms[0], molecule.atoms.len(), bond_index)?;
let at1 = source_atom_index(bond.atoms[1], molecule.atoms.len(), bond_index)?;
if at0 >= nlimit || at1 >= nlimit {
continue;
}
neighbours[at0].push(at1, bond_index, at0)?;
neighbours[at1].push(at0, bond_index, at1)?;
}
Ok(neighbours)
}
fn source_atom_index(
source_index: i32,
atom_count: usize,
bond_index: usize,
) -> Result<usize, FingerprintError> {
let index = source_index
.checked_sub(1)
.and_then(|value| usize::try_from(value).ok());
index
.filter(|&value| value < atom_count)
.ok_or_else(|| FingerprintError::AvalonConversion {
reason: format!(
"Avalon bond {} references invalid atom {}",
bond_index + 1,
source_index
),
})
}
fn implicit_hydrogens(
symbol: &str,
nsingle: i32,
naromatic: i32,
ndouble: i32,
ntriple: i32,
radical: i32,
charge: i32,
) -> i32 {
let mut bond_electrons = nsingle + 2 * ndouble + 3 * ntriple;
if radical != 0 {
bond_electrons += 1;
}
bond_electrons += match naromatic {
0 => 0,
1 => 2,
2 => 3,
3 => 4,
count => count + 1,
};
for entry in VALENCE_TABLE
.iter()
.filter(|entry| entry.atom_type == symbol)
{
if charge == 0 {
let mut value = entry.from_valence;
while value <= entry.to_valence {
let hydrogens = value - bond_electrons;
if hydrogens >= 0 {
return hydrogens;
}
value += entry.step_valence;
}
} else if charge > 0 {
let mut value = entry.from_valence;
while value <= entry.to_valence {
let hydrogens = value - bond_electrons + charge;
if hydrogens < 0 {
value += entry.step_valence;
continue;
}
return if entry.lone_pairs > 0 { hydrogens } else { 0 };
}
} else {
let mut value = entry.from_valence;
while value <= entry.to_valence {
let mut hydrogens = value - bond_electrons - charge;
if hydrogens < 0 {
value += entry.step_valence;
continue;
}
if entry.pair_deficit < hydrogens {
hydrogens = entry.pair_deficit;
}
return if entry.lone_pairs > 0 { 0 } else { hydrogens };
}
}
}
0
}
fn compute_implicit_h(
molecule: &MoleculeState,
h_count: &mut [i32],
) -> Result<(), FingerprintError> {
if h_count.len() != molecule.atoms.len() + 1 {
return Err(FingerprintError::AvalonConversion {
reason: "Avalon H-count array has the wrong one-based length".to_string(),
});
}
let array_len = molecule.atoms.len() + 1;
let mut single_bond = vec![0_i32; array_len];
let mut aromatic_bond = vec![0_i32; array_len];
let mut double_bond = vec![0_i32; array_len];
let mut triple_bond = vec![0_i32; array_len];
let mut radical = vec![0_i32; array_len];
let mut charge = vec![0_i32; array_len];
for (index, atom) in molecule.atoms.iter().enumerate() {
radical[index + 1] = i32::from(atom.radical == DOUBLET);
charge[index + 1] = atom.charge;
}
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
let at0 = source_atom_number(bond.atoms[0], molecule.atoms.len(), bond_index)?;
let at1 = source_atom_number(bond.atoms[1], molecule.atoms.len(), bond_index)?;
match bond.bond_type {
SINGLE => {
single_bond[at0] += 1;
single_bond[at1] += 1;
}
DOUBLE => {
double_bond[at0] += 1;
double_bond[at1] += 1;
}
TRIPLE => {
triple_bond[at0] += 1;
triple_bond[at1] += 1;
}
AROMATIC => {
aromatic_bond[at0] += 1;
aromatic_bond[at1] += 1;
}
_ => {
single_bond[at0] += 1;
single_bond[at1] += 1;
} }
}
for (index, atom) in molecule.atoms.iter().enumerate() {
if h_count[index + 1] == 0 {
h_count[index + 1] = implicit_hydrogens(
&atom.atom_symbol,
single_bond[index + 1],
aromatic_bond[index + 1],
double_bond[index + 1],
triple_bond[index + 1],
radical[index + 1],
charge[index + 1],
);
if h_count[index + 1] < 0 {
h_count[index + 1] = 0;
}
}
}
Ok(())
}
fn source_atom_number(
source_index: i32,
atom_count: usize,
bond_index: usize,
) -> Result<usize, FingerprintError> {
usize::try_from(source_index)
.ok()
.filter(|&value| value > 0 && value <= atom_count)
.ok_or_else(|| FingerprintError::AvalonConversion {
reason: format!(
"Avalon bond {} references invalid atom {}",
bond_index + 1,
source_index
),
})
}
fn guess_h_counts_from_substitution(
molecule: &mut MoleculeState,
neighbours: &[Neighbourhood],
) -> Result<(), FingerprintError> {
if neighbours.len() != molecule.atoms.len() {
return Err(FingerprintError::AvalonConversion {
reason: "Avalon neighbourhood array has the wrong length".to_string(),
});
}
for (index, neighbourhood) in neighbours.iter().enumerate() {
let atom = &molecule.atoms[index];
if atom.query_h_count != NONE {
continue;
}
if atom.charge != 0 {
continue;
}
if atom.radical != 0 {
continue;
}
let (mut nsingle, mut ndouble, mut naromatic, mut ntriple, mut nother) =
(0_i32, 0_i32, 0_i32, 0_i32, 0_i32);
let mut nexplicit = 0_i32;
for &neighbour_atom in neighbourhood.atoms() {
if matches!(
molecule.atoms[neighbour_atom].atom_symbol.as_str(),
"H" | "D" | "T"
) {
nexplicit += 1;
}
}
for &bond_index in neighbourhood.bonds() {
match molecule.bonds[bond_index].bond_type {
SINGLE => nsingle += 1,
DOUBLE => ndouble += 1,
AROMATIC => naromatic += 1,
TRIPLE => ntriple += 1,
_ => nother += 1,
}
}
if nother > 0 {
continue;
}
let atom = &mut molecule.atoms[index];
if atom.atom_symbol == "C" {
if naromatic == 2 {
nsingle += 1;
ndouble += 1;
naromatic = 0;
}
if naromatic > 0 {
continue;
}
if atom.sub_desc == NONE {
if nsingle + 2 * ndouble + 3 * ntriple == 4 {
atom.sub_desc = SUB_AS_IS;
}
} else if atom.sub_desc == SUB_AS_IS {
atom.query_h_count =
nexplicit + ZERO_COUNT + 4 - nsingle - 2 * ndouble - 3 * ntriple;
}
} else if atom.atom_symbol == "O" {
if ntriple + naromatic != 0 {
continue;
}
if atom.sub_desc == SUB_AS_IS {
atom.query_h_count = nexplicit + ZERO_COUNT + 2 - nsingle - 2 * ndouble;
}
} else if atom.atom_symbol == "N" {
if naromatic != 0 {
continue;
}
if atom.sub_desc == SUB_AS_IS {
atom.query_h_count = nexplicit + ZERO_COUNT + 3 - nsingle - 2 * ndouble;
}
}
}
Ok(())
}
pub(super) fn collect_hydrogen_counts(
molecule: &mut MoleculeState,
neighbours: &[Neighbourhood],
as_query: bool,
) -> Result<Vec<i32>, FingerprintError> {
let mut h_count = vec![0_i32; molecule.atoms.len() + 1];
if as_query {
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
let at0 = source_atom_number(bond.atoms[0], molecule.atoms.len(), bond_index)?;
let at1 = source_atom_number(bond.atoms[1], molecule.atoms.len(), bond_index)?;
if molecule.atoms[at0 - 1].atom_symbol == "H" {
h_count[at1] += 1;
} else if molecule.atoms[at1 - 1].atom_symbol == "H" {
h_count[at0] += 1;
}
}
guess_h_counts_from_substitution(molecule, neighbours)?;
for (index, atom) in molecule.atoms.iter().enumerate() {
if atom.query_h_count != NONE {
h_count[index + 1] = atom.query_h_count - ZERO_COUNT;
}
}
} else {
compute_implicit_h(molecule, &mut h_count)?;
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
let at0 = source_atom_number(bond.atoms[0], molecule.atoms.len(), bond_index)?;
let at1 = source_atom_number(bond.atoms[1], molecule.atoms.len(), bond_index)?;
if molecule.atoms[at0 - 1].atom_symbol == "H" {
h_count[at1] += 1;
} else if molecule.atoms[at1 - 1].atom_symbol == "H" {
h_count[at0] += 1;
}
}
}
Ok(h_count)
}
pub(super) fn ring_state(
molecule: &MoleculeState,
) -> Result<(Vec<i32>, Vec<i32>), FingerprintError> {
let mut atom_status = vec![0_i32; molecule.atoms.len()];
let mut bond_status = vec![0_i32; molecule.bonds.len()];
if molecule.bonds.is_empty() {
return Ok((atom_status, bond_status));
}
let mut bonds = Vec::with_capacity(molecule.bonds.len());
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
source_atom_number(bond.atoms[0], molecule.atoms.len(), bond_index)?;
source_atom_number(bond.atoms[1], molecule.atoms.len(), bond_index)?;
bonds.push([bond.atoms[0] as usize, bond.atoms[1] as usize]);
}
let mut rings = ring_list(&bonds);
combine_rings(&mut rings);
for ring in &rings {
for (bond_index, status) in bond_status.iter_mut().enumerate() {
if ring.bond_set.contains(bond_index) {
*status += 1;
}
}
}
for (bond_index, bond) in molecule.bonds.iter().enumerate() {
if bond_status[bond_index] > 0 {
atom_status[bond.atoms[0] as usize - 1] += 1;
atom_status[bond.atoms[1] as usize - 1] += 1;
}
}
Ok((atom_status, bond_status))
}
#[allow(clippy::too_many_arguments)]
fn mark_recursive(
molecule: &mut MoleculeState,
touched_atoms: &mut [bool],
touched_bonds: &mut [bool],
start_index: usize,
path_length: usize,
current_index: usize,
max_size: usize,
neighbours: &[Neighbourhood],
) {
for ligand in 0..neighbours[current_index].n_ligands {
let atom_index = neighbours[current_index].atoms[ligand];
if atom_index < start_index {
continue;
}
if atom_index == start_index {
if path_length < 3 {
continue;
}
for (index, &touched) in touched_atoms.iter().enumerate() {
if touched {
molecule.atoms[index].rsize_flags |= 1_u32 << path_length;
}
}
for (index, &touched) in touched_bonds.iter().enumerate() {
if touched {
molecule.bonds[index].rsize_flags |= 1_u32 << path_length;
}
}
continue;
}
if touched_atoms[atom_index] {
continue;
}
if path_length + 1 > max_size {
continue;
}
if molecule.atoms[atom_index].rsize_flags == 0 {
continue;
}
let bond_index = neighbours[current_index].bonds[ligand];
if molecule.bonds[bond_index].rsize_flags == 0 {
continue;
}
touched_atoms[atom_index] = true;
touched_bonds[bond_index] = true;
mark_recursive(
molecule,
touched_atoms,
touched_bonds,
start_index,
path_length + 1,
atom_index,
max_size,
neighbours,
);
touched_atoms[atom_index] = false;
touched_bonds[bond_index] = false;
}
}
pub(super) fn set_ring_size_flags(
molecule: &mut MoleculeState,
max_size: usize,
neighbours: &[Neighbourhood],
) -> Result<(), FingerprintError> {
if neighbours.len() != molecule.atoms.len() {
return Err(FingerprintError::AvalonConversion {
reason: "Avalon neighbourhood array has the wrong length".to_string(),
});
}
if max_size >= u32::BITS as usize {
return Err(FingerprintError::AvalonConversion {
reason: "Avalon ring-size bit limit exceeds the source field width".to_string(),
});
}
let (atom_status, bond_status) = ring_state(molecule)?;
for (atom, status) in molecule.atoms.iter_mut().zip(atom_status) {
atom.rsize_flags = u32::from(status > 0);
}
for (bond, status) in molecule.bonds.iter_mut().zip(bond_status) {
bond.rsize_flags = u32::from(status > 0);
}
let mut touched_atoms = vec![false; molecule.atoms.len()];
let mut touched_bonds = vec![false; molecule.bonds.len()];
for index in 0..molecule.atoms.len() {
if molecule.atoms[index].rsize_flags == 0 {
continue;
}
touched_atoms[index] = true;
mark_recursive(
molecule,
&mut touched_atoms,
&mut touched_bonds,
index,
1,
index,
max_size,
neighbours,
);
touched_atoms[index] = false;
}
Ok(())
}
#[cfg(test)]
mod tests {
use super::*;
use crate::properties::avalon_fingerprint::reaccs::{Atom, Bond};
fn state_with_bonds(atom_count: usize, endpoints: &[[i32; 2]]) -> MoleculeState {
MoleculeState {
atoms: vec![Atom::default(); atom_count],
bonds: endpoints
.iter()
.map(|&atoms| Bond {
atoms,
..Bond::default()
})
.collect(),
..MoleculeState::default()
}
}
fn atom(symbol: &str) -> Atom {
Atom {
atom_symbol: symbol.to_string(),
..Atom::default()
}
}
fn bond(atoms: [i32; 2], bond_type: i32) -> Bond {
Bond {
atoms,
bond_type,
..Bond::default()
}
}
#[test]
fn setup_neighbourhood_preserves_source_bond_table_order() {
let molecule = state_with_bonds(4, &[[1, 3], [2, 1], [1, 4], [3, 2]]);
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
assert_eq!(neighbours[0].atoms(), &[2, 1, 3]);
assert_eq!(neighbours[0].bonds(), &[0, 1, 2]);
assert_eq!(neighbours[1].atoms(), &[0, 2]);
assert_eq!(neighbours[1].bonds(), &[1, 3]);
assert_eq!(neighbours[2].atoms(), &[0, 1]);
assert_eq!(neighbours[2].bonds(), &[0, 3]);
assert_eq!(neighbours[3].atoms(), &[0]);
assert_eq!(neighbours[3].bonds(), &[2]);
}
#[test]
fn setup_neighbourhood_applies_source_nlimit_to_both_endpoints() {
let molecule = state_with_bonds(4, &[[1, 2], [2, 3], [1, 4]]);
let neighbours = setup_neighbourhood(&molecule, 2).unwrap();
assert_eq!(neighbours[0].atoms(), &[1]);
assert_eq!(neighbours[1].atoms(), &[0]);
assert!(neighbours[2].atoms().is_empty());
assert!(neighbours[3].atoms().is_empty());
}
#[test]
fn setup_neighbourhood_accepts_twenty_and_rejects_twenty_one_ligands() {
let mut endpoints = Vec::new();
for atom in 2..=21 {
endpoints.push([1, atom]);
}
let molecule = state_with_bonds(21, &endpoints);
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
assert_eq!(neighbours[0].atoms().len(), MAX_NEIGHBOURS);
endpoints.push([1, 22]);
let molecule = state_with_bonds(22, &endpoints);
assert!(matches!(
setup_neighbourhood(&molecule, molecule.atoms.len()),
Err(FingerprintError::AvalonConversion { reason })
if reason == "Avalon neighbourhood capacity exceeded at atom 1"
));
}
#[test]
fn setup_neighbourhood_rejects_invalid_one_based_atom_indices() {
for invalid in [0, -1, 3] {
let molecule = state_with_bonds(2, &[[1, invalid]]);
assert!(matches!(
setup_neighbourhood(&molecule, molecule.atoms.len()),
Err(FingerprintError::AvalonConversion { .. })
));
}
}
#[test]
fn implicit_hydrogens_uses_source_valence_and_aromatic_electron_rules() {
assert_eq!(implicit_hydrogens("C", 0, 0, 0, 0, 0, 0), 4);
assert_eq!(implicit_hydrogens("C", 1, 0, 0, 0, 0, 0), 3);
assert_eq!(implicit_hydrogens("N", 0, 0, 0, 0, 0, 0), 3);
assert_eq!(implicit_hydrogens("N", 0, 2, 0, 0, 0, 0), 0);
assert_eq!(implicit_hydrogens("C", 0, 4, 0, 0, 0, 0), 0);
assert_eq!(implicit_hydrogens("Xe", 0, 0, 0, 0, 0, 0), 0);
}
#[test]
fn compute_implicit_h_preserves_nonzero_one_based_slots() {
let molecule = MoleculeState {
atoms: vec![atom("C"), atom("C")],
bonds: vec![bond([1, 2], SINGLE)],
..MoleculeState::default()
};
let mut h_count = vec![77, 2, 0];
compute_implicit_h(&molecule, &mut h_count).unwrap();
assert_eq!(h_count, vec![77, 2, 3]);
}
#[test]
fn ordinary_counts_add_explicit_h_after_implicit_h() {
let mut molecule = MoleculeState {
atoms: vec![atom("C"), atom("H")],
bonds: vec![bond([1, 2], SINGLE)],
..MoleculeState::default()
};
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
let h_count = collect_hydrogen_counts(&mut molecule, &neighbours, false).unwrap();
assert_eq!(h_count, vec![0, 4, 0]);
}
#[test]
fn query_h_count_overrides_precounted_explicit_h() {
let mut query_atom = atom("N");
query_atom.query_h_count = ZERO_COUNT + 3;
let mut molecule = MoleculeState {
atoms: vec![query_atom, atom("H")],
bonds: vec![bond([1, 2], SINGLE)],
..MoleculeState::default()
};
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
let h_count = collect_hydrogen_counts(&mut molecule, &neighbours, true).unwrap();
assert_eq!(h_count, vec![0, 3, 0]);
}
#[test]
fn substitution_guess_counts_hydrogen_isotopes_but_rejects_query_bonds() {
let mut oxygen = atom("O");
oxygen.sub_desc = SUB_AS_IS;
let mut molecule = MoleculeState {
atoms: vec![oxygen, atom("C"), atom("D")],
bonds: vec![bond([1, 2], SINGLE), bond([1, 3], SINGLE)],
..MoleculeState::default()
};
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
guess_h_counts_from_substitution(&mut molecule, &neighbours).unwrap();
assert_eq!(molecule.atoms[0].query_h_count, ZERO_COUNT + 1);
molecule.atoms[0].query_h_count = NONE;
molecule.bonds[1].bond_type = 5;
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
guess_h_counts_from_substitution(&mut molecule, &neighbours).unwrap();
assert_eq!(molecule.atoms[0].query_h_count, NONE);
}
#[test]
fn carbon_none_substitution_only_transitions_to_as_is_on_first_pass() {
let mut molecule = MoleculeState {
atoms: vec![atom("C"), atom("C"), atom("C"), atom("C"), atom("C")],
bonds: vec![
bond([1, 2], SINGLE),
bond([1, 3], SINGLE),
bond([1, 4], SINGLE),
bond([1, 5], SINGLE),
],
..MoleculeState::default()
};
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
guess_h_counts_from_substitution(&mut molecule, &neighbours).unwrap();
assert_eq!(molecule.atoms[0].sub_desc, SUB_AS_IS);
assert_eq!(molecule.atoms[0].query_h_count, NONE);
}
#[test]
fn ring_state_counts_basis_membership_and_attached_ring_bonds() {
let molecule = MoleculeState {
atoms: vec![atom("C"); 4],
bonds: vec![
bond([1, 2], SINGLE),
bond([2, 3], SINGLE),
bond([3, 4], SINGLE),
bond([4, 1], SINGLE),
bond([1, 3], SINGLE),
],
..MoleculeState::default()
};
let (atom_status, bond_status) = ring_state(&molecule).unwrap();
assert_eq!(atom_status, vec![3, 2, 3, 2]);
assert_eq!(bond_status, vec![1, 1, 1, 1, 2]);
}
#[test]
fn ring_size_flags_enumerate_triangles_and_outer_square() {
let mut molecule = MoleculeState {
atoms: vec![atom("C"); 4],
bonds: vec![
bond([1, 2], SINGLE),
bond([2, 3], SINGLE),
bond([3, 4], SINGLE),
bond([4, 1], SINGLE),
bond([1, 3], SINGLE),
],
..MoleculeState::default()
};
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
set_ring_size_flags(&mut molecule, 14, &neighbours).unwrap();
let triangle_and_square = 1 | (1 << 3) | (1 << 4);
assert!(
molecule
.atoms
.iter()
.all(|atom| atom.rsize_flags == triangle_and_square)
);
assert_eq!(
molecule
.bonds
.iter()
.map(|bond| bond.rsize_flags)
.collect::<Vec<_>>(),
vec![
triangle_and_square,
triangle_and_square,
triangle_and_square,
triangle_and_square,
1 | (1 << 3),
]
);
}
#[test]
fn ring_size_flags_clear_stale_flags_on_acyclic_graph() {
let mut molecule = MoleculeState {
atoms: vec![atom("C"); 3],
bonds: vec![bond([1, 2], SINGLE), bond([2, 3], SINGLE)],
..MoleculeState::default()
};
for atom in &mut molecule.atoms {
atom.rsize_flags = u32::MAX;
}
for bond in &mut molecule.bonds {
bond.rsize_flags = u32::MAX;
}
let neighbours = setup_neighbourhood(&molecule, molecule.atoms.len()).unwrap();
set_ring_size_flags(&mut molecule, 14, &neighbours).unwrap();
assert!(molecule.atoms.iter().all(|atom| atom.rsize_flags == 0));
assert!(molecule.bonds.iter().all(|bond| bond.rsize_flags == 0));
}
}