use std::collections::VecDeque;
use crate::chemistry::valence::rdkit_atomic_mass;
use crate::{
Atom, AtomId, Bond, BondId, BondOrder, BondStereo, ChiralTag, Molecule, RingInfo,
ValenceAssignment, ValenceModel,
};
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
#[allow(non_camel_case_types)]
pub(crate) enum Descriptor {
None,
Unknown,
ns,
R,
S,
r,
s,
seqTrans,
seqCis,
E,
Z,
M,
P,
m,
p,
SP_4,
TBPY_5,
OC_6,
}
impl Descriptor {
pub(crate) const ALL_IN_RDKIT_ORDER: [Self; 18] = [
Self::None,
Self::Unknown,
Self::ns,
Self::R,
Self::S,
Self::r,
Self::s,
Self::seqTrans,
Self::seqCis,
Self::E,
Self::Z,
Self::M,
Self::P,
Self::m,
Self::p,
Self::SP_4,
Self::TBPY_5,
Self::OC_6,
];
}
pub(crate) fn descriptor_to_string(desc: Descriptor) -> &'static str {
match desc {
Descriptor::None => "NONE",
Descriptor::Unknown => "UNKNOWN",
Descriptor::ns => "ns",
Descriptor::R => "R",
Descriptor::S => "S",
Descriptor::r => "r",
Descriptor::s => "s",
Descriptor::seqTrans => "e",
Descriptor::seqCis => "z",
Descriptor::E => "E",
Descriptor::Z => "Z",
Descriptor::M => "M",
Descriptor::P => "P",
Descriptor::m => "m",
Descriptor::p => "p",
Descriptor::SP_4 => "SP_4",
Descriptor::TBPY_5 => "TBPY_5",
Descriptor::OC_6 => "OC_6",
}
}
#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
pub(crate) enum CipLabelerError {
#[error("CIPLabeler atom index {index} is out of range for {atom_count} atoms")]
AtomIndexOutOfRange { index: usize, atom_count: usize },
#[error("CIPLabeler bond index {index} is out of range for {bond_count} bonds")]
BondIndexOutOfRange { index: usize, bond_count: usize },
#[error("CIPLabeler bond {bond} is not incident to atom {atom}")]
BondNotIncident { bond: usize, atom: usize },
#[error("CIPLabeler non integer-order bond is not allowed: {order:?}")]
NonIntegerBondOrder { order: BondOrder },
#[error("CIPLabeler node {node} is not an endpoint of edge {edge}")]
EdgeEndpointMismatch { edge: usize, node: usize },
#[error("Digraph generation failed: more than {limit} nodes found.")]
TooManyNodes { limit: usize },
#[error("Max Iterations Exceeded in CIP label calculation")]
MaxIterationsExceeded,
#[error("CIPLabeler unexpected up-edge ordering")]
UnexpectedUpEdgeOrdering,
#[error("No sequence rule provided")]
NoSequenceRuleProvided,
#[error("Descriptor lists should be the same length!")]
DescriptorListLengthMismatch,
#[error("Invalid stereo descriptor")]
InvalidStereoDescriptor,
#[error("Substituents should be topologically equivalent!")]
SubstituentsShouldBeTopologicallyEquivalent,
#[error("Something unexpected!")]
SomethingUnexpected,
#[error("Rule4b instance not in rule set")]
Rule4bInstanceNotInRuleSet,
#[error("Rule5New instance not in rule set")]
Rule5NewInstanceNotInRuleSet,
#[error("Parity vectors must have size 4.")]
ParityVectorsMustHaveSize4,
#[error("CIPLabeler Configuration requires at least one focus atom")]
EmptyConfigurationFoci,
#[error("CIPLabeler Tetrahedral received a bad atom")]
BadTetrahedralAtom,
#[error("CIPLabeler Tetrahedral received a bad config")]
BadTetrahedralConfig,
#[error("CIPLabeler Tetrahedral configuration must have 4 carriers")]
TetrahedralConfigurationMustHave4Carriers,
#[error("Received a Descriptor that is not supported for atoms")]
DescriptorNotSupportedForAtoms,
#[error("Received an invalid Atom Descriptor")]
InvalidAtomDescriptor,
#[error("Could not calculate parity! Carrier mismatch")]
CarrierMismatch,
#[error("CIPLabeler Sp2Bond received bad foci")]
BadSp2BondFoci,
#[error("CIPLabeler Sp2Bond received bad config")]
BadSp2BondConfig,
#[error("CIPLabeler Sp2Bond has incorrect number of stereo atoms")]
IncorrectNumberOfStereoAtoms,
#[error("Received a Descriptor that is not supported for double bonds")]
DescriptorNotSupportedForDoubleBonds,
#[error("Received an invalid Bond Descriptor")]
InvalidBondDescriptor,
#[error("CIPLabeler AtropisomerBond received bad foci")]
BadAtropisomerBondFoci,
#[error("CIPLabeler AtropisomerBond received bad config")]
BadAtropisomerBondConfig,
#[error("Received a Descriptor that is not supported for atropisomer bonds")]
DescriptorNotSupportedForAtropisomerBonds,
#[error("CIPLabeler dependency is not ported yet: {dependency}")]
UnsupportedDependency { dependency: &'static str },
#[error(transparent)]
RingFinding(#[from] crate::RingFindingError),
#[error(transparent)]
Valence(#[from] crate::ValenceError),
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub(crate) struct RationalI32 {
numerator: i32,
denominator: i32,
}
impl RationalI32 {
fn new(numerator: i32, denominator: i32) -> Self {
debug_assert_ne!(
denominator, 0,
"boost::rational does not accept 0 as denominator"
);
let mut numerator = numerator;
let mut denominator = denominator;
if denominator < 0 {
numerator = -numerator;
denominator = -denominator;
}
let divisor = gcd_i32(numerator, denominator);
Self {
numerator: numerator / divisor,
denominator: denominator / divisor,
}
}
fn assign(&mut self, numerator: i32, denominator: i32) {
*self = Self::new(numerator, denominator);
}
#[cfg(test)]
fn tuple(self) -> (i32, i32) {
(self.numerator, self.denominator)
}
}
fn gcd_i32(a: i32, b: i32) -> i32 {
let mut a = i64::from(a).abs();
let mut b = i64::from(b).abs();
while b != 0 {
let r = a % b;
a = b;
b = r;
}
if a == 0 { 1 } else { a as i32 }
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum MancudeType {
Cv4D3,
Nv3D2,
Nv4D3Plus,
Nv2D2Minus,
Cv3D3Minus,
Ov3D2Plus,
Other,
}
pub(crate) struct CipMol<'a> {
molecule: &'a Molecule,
rings: Option<RingInfo>,
kekulized_bond_orders: Option<Vec<BondOrder>>,
fractional_atomic_numbers: Option<Vec<RationalI32>>,
valence: Option<ValenceAssignment>,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub(crate) struct CipNodeId(usize);
impl CipNodeId {
fn new(index: usize) -> Self {
Self(index)
}
fn index(self) -> usize {
self.0
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub(crate) struct CipEdgeId(usize);
impl CipEdgeId {
fn new(index: usize) -> Self {
Self(index)
}
fn index(self) -> usize {
self.0
}
}
#[derive(Debug, Clone, PartialEq)]
pub(crate) struct CipEdge {
beg: CipNodeId,
end: CipNodeId,
bond_idx: Option<usize>,
aux: Descriptor,
}
impl CipEdge {
pub(crate) fn new(beg: CipNodeId, end: CipNodeId, bond_idx: Option<usize>) -> Self {
Self {
beg,
end,
bond_idx,
aux: Descriptor::None,
}
}
pub(crate) fn get_other(
&self,
self_id: CipEdgeId,
node: CipNodeId,
) -> Result<CipNodeId, CipLabelerError> {
if self.is_beg(node) {
Ok(self.get_end())
} else if self.is_end(node) {
Ok(self.get_beg())
} else {
Err(CipLabelerError::EdgeEndpointMismatch {
edge: self_id.index(),
node: node.index(),
})
}
}
pub(crate) fn get_beg(&self) -> CipNodeId {
self.beg
}
pub(crate) fn get_end(&self) -> CipNodeId {
self.end
}
pub(crate) fn get_bond_idx(&self) -> Option<usize> {
self.bond_idx
}
pub(crate) fn get_aux(&self) -> Descriptor {
self.aux
}
pub(crate) fn is_beg(&self, node: CipNodeId) -> bool {
node == self.beg
}
pub(crate) fn is_end(&self, node: CipNodeId) -> bool {
node == self.end
}
pub(crate) fn set_aux(&mut self, aux: Descriptor) {
self.aux = aux;
}
pub(crate) fn flip(&mut self) {
std::mem::swap(&mut self.beg, &mut self.end);
}
}
#[derive(Debug, Clone, PartialEq)]
pub(crate) struct CipNode {
digraph: usize,
atom_idx: Option<usize>,
distance: i32,
atomic_num_fraction: RationalI32,
atomic_mass: f64,
aux: Descriptor,
flags: i32,
edges: Vec<CipEdgeId>,
visit: Vec<i8>,
}
pub(crate) struct CipDigraph<'a> {
mol: CipMol<'a>,
origin: CipNodeId,
root: CipNodeId,
rule6_ref: Option<usize>,
atropisomer_mode: bool,
nodes: Vec<CipNode>,
edges: Vec<CipEdge>,
}
pub(crate) struct CipLabelerContext {
remaining_call_count: u32,
}
impl CipLabelerContext {
pub(crate) const CONSTITUTIONAL_RULE_TIMEOUT: u32 = 2_000;
pub(crate) fn new(max_recursive_iterations: u32) -> Self {
Self {
remaining_call_count: if max_recursive_iterations == 0 {
u32::MAX
} else {
max_recursive_iterations
},
}
}
fn with_remaining_call_count(remaining_call_count: u32) -> Self {
Self {
remaining_call_count,
}
}
fn decrement_remaining_call_count_and_check(&mut self) -> bool {
self.remaining_call_count = self.remaining_call_count.wrapping_sub(1);
self.remaining_call_count > 0
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) struct CipPriority {
unique: bool,
pseudo_asym: bool,
}
impl CipPriority {
pub(crate) fn new(unique: bool, pseudo_asym: bool) -> Self {
Self {
unique,
pseudo_asym,
}
}
pub(crate) fn is_unique(self) -> bool {
self.unique
}
pub(crate) fn is_pseudo_asymetric(self) -> bool {
self.pseudo_asym
}
}
pub(crate) trait CipSequenceRule {
fn get_bond_label(&self, edge: &CipEdge) -> Descriptor {
if edge.get_bond_idx().is_none() {
return Descriptor::None;
}
edge.get_aux()
}
fn get_comparison(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
deep: bool,
) -> Result<i32, CipLabelerError> {
self.get_comparison_with_sort_rules(None, digraph, context, a, b, deep)
}
fn get_comparison_with_sort_rules(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
deep: bool,
) -> Result<i32, CipLabelerError> {
if deep {
recursive_compare_sequence_rule(self, sort_rules, digraph, context, a, b)
} else {
self.compare_with_sort_rules(sort_rules, digraph, context, a, b)
}
}
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError>;
fn compare_with_sort_rules(
&self,
_sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
self.compare(digraph, context, a, b)
}
fn recursive_compare(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
recursive_compare_sequence_rule(self, None, digraph, context, a, b)
}
fn recursive_compare_with_sort_rules(
&self,
sort_rules: &[&dyn CipSequenceRule],
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
recursive_compare_sequence_rule(self, Some(sort_rules), digraph, context, a, b)
}
fn sort(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
edges: &mut [CipEdgeId],
deep: bool,
) -> Result<CipPriority, CipLabelerError> {
prioritize_single_sequence_rule(self, digraph, context, node, edges, deep)
}
fn are_up_edges(
&self,
digraph: &CipDigraph<'_>,
a_node: CipNodeId,
b_node: CipNodeId,
a_edge: CipEdgeId,
b_edge: CipEdgeId,
) -> Result<bool, CipLabelerError> {
if digraph.edge(a_edge).is_end(a_node) {
if !digraph.edge(b_edge).is_end(b_node) {
return Err(CipLabelerError::UnexpectedUpEdgeOrdering);
}
return Ok(true);
}
Ok(false)
}
}
fn recursive_compare_sequence_rule<R: CipSequenceRule + ?Sized>(
rule: &R,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
mut a: CipEdgeId,
mut b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
if !context.decrement_remaining_call_count_and_check() {
return Err(CipLabelerError::MaxIterationsExceeded);
}
let mut cmp = rule.compare(digraph, context, a, b)?;
if cmp != 0 {
return Ok(cmp);
}
let mut a_queue = vec![a];
let mut b_queue = vec![b];
let mut pos = 0_usize;
while pos < a_queue.len() && pos < b_queue.len() {
a = a_queue[pos];
b = b_queue[pos];
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let mut as_edges = digraph.node_edges(a_end)?;
let mut bs_edges = digraph.node_edges(b_end)?;
sort_edges_for_sequence_rule(
rule,
sort_rules,
digraph,
context,
a_end,
&mut as_edges,
false,
)?;
sort_edges_for_sequence_rule(
rule,
sort_rules,
digraph,
context,
b_end,
&mut bs_edges,
false,
)?;
let sizediff = three_way_comparison_i32(as_edges.len() as i32, bs_edges.len() as i32);
for (a_edge, b_edge) in as_edges.iter().zip(bs_edges.iter()) {
if rule.are_up_edges(digraph, a_end, b_end, *a_edge, *b_edge)? {
continue;
}
cmp = rule.compare(digraph, context, *a_edge, *b_edge)?;
if cmp != 0 {
return Ok(cmp);
}
}
if sizediff != 0 {
return Ok(sizediff);
}
sort_edges_for_sequence_rule(
rule,
sort_rules,
digraph,
context,
a_end,
&mut as_edges,
true,
)?;
sort_edges_for_sequence_rule(
rule,
sort_rules,
digraph,
context,
b_end,
&mut bs_edges,
true,
)?;
for (a_edge, b_edge) in as_edges.iter().zip(bs_edges.iter()) {
if rule.are_up_edges(digraph, a_end, b_end, *a_edge, *b_edge)? {
continue;
}
cmp = rule.compare(digraph, context, *a_edge, *b_edge)?;
if cmp != 0 {
return Ok(cmp);
}
a_queue.push(*a_edge);
b_queue.push(*b_edge);
}
pos += 1;
}
Ok(0)
}
fn sort_edges_for_sequence_rule<R: CipSequenceRule + ?Sized>(
rule: &R,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
edges: &mut [CipEdgeId],
deep: bool,
) -> Result<CipPriority, CipLabelerError> {
if let Some(sort_rules) = sort_rules {
CipSort::from_rules(sort_rules.to_vec()).prioritize(digraph, context, node, edges, deep)
} else {
rule.sort(digraph, context, node, edges, deep)
}
}
fn prioritize_single_sequence_rule<R: CipSequenceRule + ?Sized>(
rule: &R,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
edges: &mut [CipEdgeId],
deep: bool,
) -> Result<CipPriority, CipLabelerError> {
let mut unique = true;
let mut num_pseudo_asym = 0_i32;
for i in 0..edges.len() {
let mut j = i;
while j > 0 {
let cmp = compare_substituents_with_rule(
rule,
digraph,
context,
node,
edges[j - 1],
edges[j],
deep,
)?;
if !(-1..=1).contains(&cmp) {
num_pseudo_asym += 1;
}
if cmp < 0 {
edges.swap(j, j - 1);
} else {
if cmp == 0 {
unique = false;
}
break;
}
j -= 1;
}
}
Ok(CipPriority::new(unique, num_pseudo_asym == 1))
}
fn compare_substituents_with_rule<R: CipSequenceRule + ?Sized>(
rule: &R,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
a: CipEdgeId,
b: CipEdgeId,
deep: bool,
) -> Result<i32, CipLabelerError> {
let a_is_beg = digraph.edge(a).is_beg(node);
let b_is_beg = digraph.edge(b).is_beg(node);
if !a_is_beg && b_is_beg {
return Ok(1);
} else if a_is_beg && !b_is_beg {
return Ok(-1);
}
rule.get_comparison(digraph, context, a, b, deep)
}
fn three_way_comparison_i32(x: i32, y: i32) -> i32 {
if x < y {
-1
} else if x == y {
0
} else {
1
}
}
pub(crate) struct CipSort<'r> {
rules: Vec<&'r dyn CipSequenceRule>,
}
impl<'r> CipSort<'r> {
pub(crate) fn new(rule: &'r dyn CipSequenceRule) -> Self {
Self { rules: vec![rule] }
}
pub(crate) fn from_rules(rules: Vec<&'r dyn CipSequenceRule>) -> Self {
Self { rules }
}
pub(crate) fn get_rules(&self) -> &[&'r dyn CipSequenceRule] {
&self.rules
}
pub(crate) fn prioritize(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
edges: &mut [CipEdgeId],
deep: bool,
) -> Result<CipPriority, CipLabelerError> {
let mut unique = true;
let mut num_pseudo_asym = 0_i32;
for i in 0..edges.len() {
let mut j = i;
while j > 0 {
let cmp = self.compare_substituents(
digraph,
context,
node,
edges[j - 1],
edges[j],
deep,
)?;
if !(-1..=1).contains(&cmp) {
num_pseudo_asym += 1;
}
if cmp < 0 {
edges.swap(j, j - 1);
} else {
if cmp == 0 {
unique = false;
}
break;
}
j -= 1;
}
}
Ok(CipPriority::new(unique, num_pseudo_asym == 1))
}
pub(crate) fn get_groups(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
sorted: &[CipEdgeId],
) -> Result<Vec<Vec<CipEdgeId>>, CipLabelerError> {
let mut groups = Vec::<Vec<CipEdgeId>>::new();
let mut prev = None;
for edge in sorted {
if prev.is_none()
|| self.compare_substituents(
digraph,
context,
digraph.edge(prev.expect("checked")).get_beg(),
prev.expect("checked"),
*edge,
true,
)? != 0
{
groups.push(Vec::new());
}
prev = Some(*edge);
groups
.last_mut()
.expect("RDKit Sort::getGroups creates a group before push")
.push(*edge);
}
Ok(groups)
}
fn compare_substituents(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
a: CipEdgeId,
b: CipEdgeId,
deep: bool,
) -> Result<i32, CipLabelerError> {
let a_is_beg = digraph.edge(a).is_beg(node);
let b_is_beg = digraph.edge(b).is_beg(node);
if !a_is_beg && b_is_beg {
return Ok(1);
} else if a_is_beg && !b_is_beg {
return Ok(-1);
}
for rule in &self.rules {
let cmp = rule.get_comparison_with_sort_rules(
Some(&self.rules),
digraph,
context,
a,
b,
deep,
)?;
if cmp != 0 {
return Ok(cmp);
}
}
Ok(0)
}
}
pub(crate) struct CipRules {
rules: Vec<Box<dyn CipSequenceRule>>,
}
impl CipRules {
pub(crate) fn new(rules: Vec<Box<dyn CipSequenceRule>>) -> Result<Self, CipLabelerError> {
let mut result = Self { rules: Vec::new() };
for rule in rules {
result.add(Some(rule))?;
}
Ok(result)
}
fn add(&mut self, rule: Option<Box<dyn CipSequenceRule>>) -> Result<(), CipLabelerError> {
let rule = rule.ok_or(CipLabelerError::NoSequenceRuleProvided)?;
self.rules.push(rule);
Ok(())
}
pub(crate) fn get_num_sub_rules(&self) -> usize {
self.rules.len()
}
pub(crate) fn get_sorter(&self) -> CipSort<'_> {
CipSort::new(self)
}
fn rule_refs(&self) -> Vec<&dyn CipSequenceRule> {
self.rules
.iter()
.map(|rule| rule.as_ref() as &dyn CipSequenceRule)
.collect()
}
}
fn cip_all_rules() -> Result<CipRules, CipLabelerError> {
CipRules::new(vec![
Box::new(CipRule1a),
Box::new(CipRule1b),
Box::new(CipRule2),
Box::new(CipRule3),
Box::new(CipRule4a),
Box::new(CipRule4b::new()),
Box::new(CipRule4c),
Box::new(CipRule5New::new()),
Box::new(CipRule6),
])
}
fn cip_constitutional_rules() -> Result<CipRules, CipLabelerError> {
CipRules::new(vec![
Box::new(CipRule1a),
Box::new(CipRule1b),
Box::new(CipRule2),
])
}
fn cip_find_configs<'a>(
molecule: &'a Molecule,
atom_mask: &[bool],
bond_mask: &[bool],
) -> Result<Vec<CipConfig<'a>>, CipLabelerError> {
let mut configs = Vec::new();
for (index, selected) in atom_mask.iter().copied().enumerate() {
if !selected {
continue;
}
let atom = molecule
.atoms()
.get(index)
.ok_or(CipLabelerError::AtomIndexOutOfRange {
index,
atom_count: molecule.num_atoms(),
})?;
if matches!(
atom.chiral_tag(),
ChiralTag::TetrahedralCw | ChiralTag::TetrahedralCcw
) {
configs.push(CipConfig::Tetrahedral(CipTetrahedral::new(
molecule, index,
)?));
}
}
for (index, selected) in bond_mask.iter().copied().enumerate() {
if !selected {
continue;
}
let bond = molecule
.bonds()
.get(index)
.ok_or(CipLabelerError::BondIndexOutOfRange {
index,
bond_count: molecule.num_bonds(),
})?;
let bond_cfg = match bond.stereo() {
BondStereo::E => BondStereo::Trans,
BondStereo::Z => BondStereo::Cis,
other => other,
};
match bond_cfg {
BondStereo::Trans | BondStereo::Cis => {
configs.push(CipConfig::Sp2Bond(CipSp2Bond::new(
molecule,
index,
bond.begin().index(),
bond.end().index(),
bond_cfg,
)?));
}
BondStereo::AtropCcw | BondStereo::AtropCw => {
configs.push(CipConfig::AtropisomerBond(CipAtropisomerBond::new(
molecule,
index,
bond.begin().index(),
bond.end().index(),
bond_cfg,
)?));
}
_ => {}
}
}
Ok(configs)
}
fn cip_label_with_center_digraph(
configs: &mut [CipConfig<'_>],
center_idx: usize,
target_idx: usize,
node: CipNodeId,
rules: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
debug_assert_ne!(center_idx, target_idx);
if center_idx < target_idx {
let (left, right) = configs.split_at_mut(target_idx);
let center = &mut left[center_idx];
let target = &mut right[0];
let digraph = center.get_digraph_mut();
target.label_with_external_digraph(node, digraph, rules, context)
} else {
let (left, right) = configs.split_at_mut(center_idx);
let target = &mut left[target_idx];
let center = &mut right[0];
let digraph = center.get_digraph_mut();
target.label_with_external_digraph(node, digraph, rules, context)
}
}
fn cip_set_center_node_aux(
configs: &mut [CipConfig<'_>],
center_idx: usize,
node: CipNodeId,
desc: Descriptor,
) {
let center = &mut configs[center_idx];
center.get_digraph_mut().nodes[node.index()].set_aux(desc);
}
fn cip_label_aux(
configs: &mut [CipConfig<'_>],
rules: &CipRules,
center_idx: usize,
context: &mut CipLabelerContext,
) -> Result<bool, CipLabelerError> {
let config_foci = configs
.iter()
.enumerate()
.filter(|(idx, _)| *idx != center_idx)
.map(|(idx, config)| (idx, config.get_foci().to_vec()))
.collect::<Vec<_>>();
let mut aux = Vec::<(CipNodeId, usize, i32)>::new();
{
let digraph = configs[center_idx].get_digraph_mut();
for (config_idx, foci) in config_foci {
for node in digraph.get_nodes(foci[0])? {
if digraph.node(node).is_duplicate() {
continue;
}
let mut low = node;
if foci.len() == 2 {
for edge in digraph.node_edges_for_atom(node, Some(foci[1]))? {
let other_node = digraph.edge(edge).get_other(edge, node)?;
if digraph.node(other_node).get_distance()
< digraph.node(node).get_distance()
{
low = other_node;
}
}
}
if !digraph.node(low).is_duplicate() {
aux.push((low, config_idx, digraph.node(low).get_distance()));
}
}
}
}
aux.sort_by(|left, right| right.2.cmp(&left.2));
let mut queue = Vec::<(CipNodeId, Descriptor)>::new();
let mut prev = i32::MAX;
for (node, config_idx, distance) in aux {
if distance < prev {
for (queued_node, desc) in queue.drain(..) {
cip_set_center_node_aux(configs, center_idx, queued_node, desc);
}
prev = distance;
}
let label =
cip_label_with_center_digraph(configs, center_idx, config_idx, node, rules, context)?;
if !queue.iter().any(|(queued_node, _)| *queued_node == node) {
queue.push((node, label));
}
}
for (queued_node, desc) in queue {
cip_set_center_node_aux(configs, center_idx, queued_node, desc);
}
Ok(true)
}
fn cip_label(
configs: &mut [CipConfig<'_>],
max_recursive_iterations: u32,
) -> Result<(), CipLabelerError> {
let constitutional_rules = cip_constitutional_rules()?;
for conf in configs.iter_mut() {
conf.reset_primary_label();
let mut context = CipLabelerContext::with_remaining_call_count(
CipLabelerContext::CONSTITUTIONAL_RULE_TIMEOUT,
);
match conf.label(&constitutional_rules, &mut context) {
Ok(desc) if desc != Descriptor::Unknown => conf.set_primary_label(desc)?,
Ok(_) | Err(CipLabelerError::MaxIterationsExceeded) => {}
Err(err) => return Err(err),
}
}
let constitutional_rules = cip_constitutional_rules()?;
let all_rules = cip_all_rules()?;
let mut context = CipLabelerContext::new(max_recursive_iterations);
for idx in 0..configs.len() {
if configs[idx].has_primary_label() {
continue;
}
let desc = configs[idx].label(&constitutional_rules, &mut context)?;
if desc != Descriptor::Unknown {
configs[idx].set_primary_label(desc)?;
} else if cip_label_aux(configs, &all_rules, idx, &mut context)? {
let desc = configs[idx].label(&all_rules, &mut context)?;
if desc != Descriptor::Unknown {
configs[idx].set_primary_label(desc)?;
}
}
}
Ok(())
}
fn cip_neighbor_order_value(indices: &[usize]) -> String {
let body = indices
.iter()
.map(usize::to_string)
.collect::<Vec<_>>()
.join(",");
format!("[{body}]")
}
fn cip_clear_selected_labels(molecule: &mut Molecule, atom_mask: &[bool], bond_mask: &[bool]) {
for (idx, selected) in atom_mask.iter().copied().enumerate() {
if !selected {
continue;
}
if let Some(atom) = molecule.topology_block_mut().atoms.get_mut(idx) {
if matches!(
atom.chiral_tag(),
ChiralTag::TetrahedralCw | ChiralTag::TetrahedralCcw
) {
atom.clear_prop("_CIPCode");
atom.clear_prop("_CIPNeighborOrder");
}
}
}
for (idx, selected) in bond_mask.iter().copied().enumerate() {
if !selected {
continue;
}
if let Some(bond) = molecule.topology_block_mut().bonds.get_mut(idx) {
if matches!(
bond.stereo(),
BondStereo::E
| BondStereo::Z
| BondStereo::Cis
| BondStereo::Trans
| BondStereo::AtropCcw
| BondStereo::AtropCw
) {
bond.clear_prop("_CIPCode");
bond.clear_prop("_CIPNeighborOrder");
}
}
}
}
fn cip_apply_primary_labels(
molecule: &mut Molecule,
labels: Vec<CipPrimaryLabel>,
) -> Result<(), CipLabelerError> {
let topology = molecule.topology_block_mut();
for label in labels {
match label {
CipPrimaryLabel::Atom(label) => {
let atom_count = topology.atoms.len();
let atom = topology.atoms.get_mut(label.atom_idx).ok_or(
CipLabelerError::AtomIndexOutOfRange {
index: label.atom_idx,
atom_count,
},
)?;
atom.set_prop("_CIPCode", label.cip_code);
atom.set_prop(
"_CIPNeighborOrder",
cip_neighbor_order_value(&label.cip_neighbor_order),
);
}
CipPrimaryLabel::Bond(label) => {
let bond_count = topology.bonds.len();
let bond = topology.bonds.get_mut(label.bond_idx).ok_or(
CipLabelerError::BondIndexOutOfRange {
index: label.bond_idx,
bond_count,
},
)?;
bond.set_stereo_atoms(Some([
AtomId::new(label.stereo_atoms[0]),
AtomId::new(label.stereo_atoms[1]),
]));
bond.set_stereo(label.stereo);
bond.set_prop("_CIPCode", label.cip_code);
bond.set_prop(
"_CIPNeighborOrder",
cip_neighbor_order_value(&label.cip_neighbor_order),
);
}
CipPrimaryLabel::AtropisomerBond(label) => {
let bond_count = topology.bonds.len();
let bond = topology.bonds.get_mut(label.bond_idx).ok_or(
CipLabelerError::BondIndexOutOfRange {
index: label.bond_idx,
bond_count,
},
)?;
bond.set_prop("_CIPCode", label.cip_code);
bond.set_prop(
"_CIPNeighborOrder",
cip_neighbor_order_value(&label.cip_neighbor_order),
);
}
}
}
Ok(())
}
pub(crate) fn assign_cip_labels_for_indices(
molecule: &Molecule,
atom_mask: &[bool],
bond_mask: &[bool],
max_recursive_iterations: u32,
) -> Result<Molecule, CipLabelerError> {
let mut labeled = molecule.clone();
labeled.properties_mut().clear_prop("_CIPComputed");
cip_clear_selected_labels(&mut labeled, atom_mask, bond_mask);
let mut configs = cip_find_configs(&labeled, atom_mask, bond_mask)?;
cip_label(&mut configs, max_recursive_iterations)?;
let labels = configs
.iter()
.filter_map(CipConfig::primary_label)
.collect::<Vec<_>>();
drop(configs);
cip_apply_primary_labels(&mut labeled, labels)?;
labeled.properties_mut().set_prop("_CIPComputed", "1");
Ok(labeled)
}
pub(crate) fn assign_cip_labels(
molecule: &Molecule,
max_recursive_iterations: u32,
) -> Result<Molecule, CipLabelerError> {
let atom_mask = vec![true; molecule.num_atoms()];
let bond_mask = vec![true; molecule.num_bonds()];
assign_cip_labels_for_indices(molecule, &atom_mask, &bond_mask, max_recursive_iterations)
}
impl CipSequenceRule for CipRules {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let sort_rules = self.rule_refs();
for rule in &sort_rules {
let value =
rule.recursive_compare_with_sort_rules(&sort_rules, digraph, context, a, b)?;
if value != 0 {
return Ok(value);
}
}
Ok(0)
}
fn get_comparison(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
_deep: bool,
) -> Result<i32, CipLabelerError> {
let sort_rules = self.rule_refs();
for rule in &sort_rules {
let value =
rule.recursive_compare_with_sort_rules(&sort_rules, digraph, context, a, b)?;
if value != 0 {
return Ok(value);
}
}
Ok(0)
}
fn sort(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
edges: &mut [CipEdgeId],
deep: bool,
) -> Result<CipPriority, CipLabelerError> {
self.get_sorter()
.prioritize(digraph, context, node, edges, deep)
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule1a;
impl CipSequenceRule for CipRule1a {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let afrac = digraph
.node(digraph.edge(a).get_end())
.get_atomic_num_fraction();
let bfrac = digraph
.node(digraph.edge(b).get_end())
.get_atomic_num_fraction();
Ok(match afrac.cmp(&bfrac) {
std::cmp::Ordering::Less => -1,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => 1,
})
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule1b;
impl CipRule1b {
const IUPAC_2013: bool = false;
}
impl CipSequenceRule for CipRule1b {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let a_node = digraph.node(a_end);
let b_node = digraph.node(b_end);
if Self::IUPAC_2013 {
return Ok(-three_way_comparison_i32(
a_node.get_distance(),
b_node.get_distance(),
));
}
if a_node.is_set(CipNode::RING_DUPLICATE) && b_node.is_set(CipNode::RING_DUPLICATE) {
return Ok(-three_way_comparison_i32(
a_node.get_distance(),
b_node.get_distance(),
));
}
if a_node.is_set(CipNode::RING_DUPLICATE) && !b_node.is_set(CipNode::RING_DUPLICATE) {
return Ok(1);
}
if !a_node.is_set(CipNode::RING_DUPLICATE) && b_node.is_set(CipNode::RING_DUPLICATE) {
return Ok(-1);
}
Ok(0)
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule2;
impl CipSequenceRule for CipRule2 {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let a_node = digraph.node(a_end);
let b_node = digraph.node(b_end);
let a_atom_num = a_node.get_atomic_num(digraph.mol())?;
let b_atom_num = b_node.get_atomic_num(digraph.mol())?;
if a_atom_num == 0 && b_atom_num == 0 {
return Ok(0);
} else if a_atom_num == 0 || b_atom_num == 0 {
return Ok(three_way_comparison_i32(
i32::from(a_atom_num),
i32::from(b_atom_num),
));
}
let a_mass_num = a_node.get_mass_num(digraph.mol())?;
let b_mass_num = b_node.get_mass_num(digraph.mol())?;
if a_mass_num == 0 && b_mass_num == 0 {
return Ok(0);
}
Ok(
match a_node
.get_atomic_mass()
.partial_cmp(&b_node.get_atomic_mass())
.expect("RDKit CIP atomic masses are finite")
{
std::cmp::Ordering::Less => -1,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => 1,
},
)
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule3;
impl CipRule3 {
fn ord(lab: Descriptor) -> i32 {
match lab {
Descriptor::E => 1,
Descriptor::Z => 2,
_ => 0,
}
}
}
impl CipSequenceRule for CipRule3 {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_ord = Self::ord(digraph.node(digraph.edge(a).get_end()).get_aux());
let b_ord = Self::ord(digraph.node(digraph.edge(b).get_end()).get_aux());
Ok(three_way_comparison_i32(a_ord, b_ord))
}
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) struct CipPairList {
descriptors: Vec<Descriptor>,
pairing: u64,
}
impl Default for CipPairList {
fn default() -> Self {
Self::new()
}
}
impl PartialOrd for CipPairList {
fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
Some(self.cmp(other))
}
}
impl Ord for CipPairList {
fn cmp(&self, other: &Self) -> std::cmp::Ordering {
match self.compare_to(other) {
Ok(-1) => std::cmp::Ordering::Less,
Ok(0) => std::cmp::Ordering::Equal,
Ok(1) => std::cmp::Ordering::Greater,
Ok(value) if value < 0 => std::cmp::Ordering::Less,
Ok(_) => std::cmp::Ordering::Greater,
Err(_) => std::cmp::Ordering::Equal,
}
}
}
impl CipPairList {
const NUM_PAIRING_BITS: usize = 64;
pub(crate) fn ref_descriptor(descriptor: Descriptor) -> Descriptor {
match descriptor {
Descriptor::R | Descriptor::M | Descriptor::seqCis => Descriptor::R,
Descriptor::S | Descriptor::P | Descriptor::seqTrans => Descriptor::S,
_ => Descriptor::None,
}
}
pub(crate) fn new() -> Self {
Self {
descriptors: Vec::new(),
pairing: 0,
}
}
pub(crate) fn with_ref(ref_descriptor: Descriptor) -> Self {
let mut result = Self::new();
result.add(ref_descriptor);
result
}
pub(crate) fn from_head_tail(head: &Self, tail: &Self) -> Self {
let mut result = Self::new();
result.add_all(&head.descriptors);
result.add_all(&tail.descriptors);
result
}
pub(crate) fn get_ref_descriptor(&self) -> Descriptor {
Self::ref_descriptor(self.descriptors[0])
}
pub(crate) fn add(&mut self, descriptor: Descriptor) -> bool {
match descriptor {
Descriptor::R
| Descriptor::S
| Descriptor::M
| Descriptor::P
| Descriptor::seqTrans
| Descriptor::seqCis => {
self.add_and_pair(descriptor);
true
}
_ => false,
}
}
pub(crate) fn add_all(&mut self, descriptors: &[Descriptor]) {
for descriptor in descriptors {
self.add(*descriptor);
}
}
pub(crate) fn get_pairing(&self) -> u64 {
self.pairing
}
pub(crate) fn compare_to(&self, that: &Self) -> Result<i32, CipLabelerError> {
if self.descriptors.len() != that.descriptors.len() {
return Err(CipLabelerError::DescriptorListLengthMismatch);
}
let this_ref = self.descriptors[0];
let that_ref = that.descriptors[0];
for i in 1..self.descriptors.len() {
if this_ref == self.descriptors[i] && that_ref != that.descriptors[i] {
return Ok(1);
}
if this_ref != self.descriptors[i] && that_ref == that.descriptors[i] {
return Ok(-1);
}
}
Ok(0)
}
pub(crate) fn to_rdkit_string(&self) -> String {
if self.descriptors.is_empty() || self.descriptors[0] == Descriptor::None {
return String::new();
}
let mut result = String::new();
let mut basis = self.descriptors[0];
result.push_str(descriptor_to_string(basis));
result.push(':');
basis = Self::ref_descriptor(basis);
for descriptor in self.descriptors.iter().skip(1) {
result.push(if basis == Self::ref_descriptor(*descriptor) {
'l'
} else {
'u'
});
}
result
}
fn add_and_pair(&mut self, descriptor: Descriptor) {
if !self.descriptors.is_empty() && self.descriptors[0] == descriptor {
self.pairing |= 1_u64 << (Self::NUM_PAIRING_BITS - 1 - self.descriptors.len());
}
self.descriptors.push(Self::ref_descriptor(descriptor));
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule4a;
impl CipRule4a {
fn ord(lab: Descriptor) -> Result<i32, CipLabelerError> {
match lab {
Descriptor::Unknown | Descriptor::ns | Descriptor::None => Ok(0),
Descriptor::r
| Descriptor::s
| Descriptor::m
| Descriptor::p
| Descriptor::E
| Descriptor::Z => Ok(1),
Descriptor::R
| Descriptor::S
| Descriptor::M
| Descriptor::P
| Descriptor::seqTrans
| Descriptor::seqCis => Ok(2),
_ => Err(CipLabelerError::InvalidStereoDescriptor),
}
}
}
impl CipSequenceRule for CipRule4a {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_ordinal = Self::ord(self.get_bond_label(digraph.edge(a)))?;
let b_ordinal = Self::ord(self.get_bond_label(digraph.edge(b)))?;
let cmp = three_way_comparison_i32(a_ordinal, b_ordinal);
if cmp != 0 {
return Ok(cmp);
}
let a_ordinal = Self::ord(digraph.node(digraph.edge(a).get_end()).get_aux())?;
let b_ordinal = Self::ord(digraph.node(digraph.edge(b).get_end()).get_aux())?;
Ok(three_way_comparison_i32(a_ordinal, b_ordinal))
}
}
#[derive(Debug, Clone, Copy)]
pub(crate) struct CipRule4b {
ref_descriptor: Descriptor,
}
impl CipRule4b {
pub(crate) fn new() -> Self {
Self {
ref_descriptor: Descriptor::None,
}
}
pub(crate) fn with_ref(ref_descriptor: Descriptor) -> Self {
Self { ref_descriptor }
}
fn get_reference_descriptors(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
node: CipNodeId,
) -> Result<Vec<Descriptor>, CipLabelerError> {
let mut result = Vec::new();
let mut prev = self.initial_level(node);
while !prev.is_empty() {
for nodes in &prev {
if self.get_reference(digraph, nodes, &mut result) {
return Ok(result);
}
}
prev = self.get_next_level(sort_rules, digraph, context, &prev)?;
}
Ok(Vec::new())
}
fn has_descriptors(
&self,
digraph: &mut CipDigraph<'_>,
node: CipNodeId,
) -> Result<bool, CipLabelerError> {
let mut queue = vec![node];
let mut pos = 0_usize;
while pos < queue.len() {
let current = queue[pos];
if digraph.node(current).get_aux() != Descriptor::None {
return Ok(true);
}
for edge_id in digraph.node_edges(current)? {
let edge = digraph.edge(edge_id);
if edge.get_end() == current {
continue;
}
if self.get_bond_label(edge) != Descriptor::None {
return Ok(true);
}
queue.push(edge.get_end());
}
pos += 1;
}
Ok(false)
}
fn get_reference(
&self,
digraph: &CipDigraph<'_>,
nodes: &[CipNodeId],
result: &mut Vec<Descriptor>,
) -> bool {
let mut right = 0_i32;
let mut left = 0_i32;
for node in nodes {
match digraph.node(*node).get_aux() {
Descriptor::None => continue,
Descriptor::R | Descriptor::M | Descriptor::seqCis => right += 1,
Descriptor::S | Descriptor::P | Descriptor::seqTrans => left += 1,
_ => {}
}
}
if right + left == 0 {
false
} else if right > left {
result.push(Descriptor::R);
true
} else if right < left {
result.push(Descriptor::S);
true
} else {
result.push(Descriptor::R);
result.push(Descriptor::S);
true
}
}
fn initial_level(&self, node: CipNodeId) -> Vec<Vec<CipNodeId>> {
vec![vec![node]]
}
fn get_next_level(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
prev_level: &[Vec<CipNodeId>],
) -> Result<Vec<Vec<CipNodeId>>, CipLabelerError> {
let mut next_level = Vec::with_capacity(4 * prev_level.len());
for prev in prev_level {
let mut tmp = Vec::<Vec<Vec<CipEdgeId>>>::new();
for node in prev {
let mut edges = digraph.non_terminal_out_edges(*node)?;
sort_edges_for_sequence_rule(
self, sort_rules, digraph, context, *node, &mut edges, true,
)?;
let groups = if let Some(sort_rules) = sort_rules {
CipSort::from_rules(sort_rules.to_vec()).get_groups(digraph, context, &edges)?
} else {
CipSort::new(self).get_groups(digraph, context, &edges)?
};
tmp.push(groups);
}
let mut size = -1_i32;
for _ in 0..tmp.len() {
let local_size = tmp.first().map(Vec::len).unwrap_or(0) as i32;
if size < 0 {
size = local_size;
} else if size != local_size {
return Err(CipLabelerError::SomethingUnexpected);
}
}
for i in 0..usize::try_from(size.max(0)).expect("nonnegative") {
let mut eq = Vec::new();
for a_tmp in &tmp {
let tmp_nodes = self.to_node_list(digraph, &a_tmp[i]);
eq.extend(tmp_nodes);
}
if !eq.is_empty() {
next_level.push(eq);
}
}
}
Ok(next_level)
}
fn to_node_list(&self, digraph: &CipDigraph<'_>, eq_edges: &[CipEdgeId]) -> Vec<CipNodeId> {
let mut eq_nodes = Vec::with_capacity(eq_edges.len());
for edge in eq_edges {
eq_nodes.push(digraph.edge(*edge).get_end());
}
eq_nodes
}
fn new_pair_lists(&self, descriptors: &[Descriptor]) -> Vec<CipPairList> {
let mut pairs = Vec::with_capacity(descriptors.len());
for descriptor in descriptors {
pairs.push(CipPairList::with_ref(*descriptor));
}
pairs
}
fn fill_pairs(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
beg: CipNodeId,
plist: &mut CipPairList,
) -> Result<(), CipLabelerError> {
let replacement_rule = CipRule4b::with_ref(plist.get_ref_descriptor());
let ref_sort_rules = self.get_ref_sorter(sort_rules, &replacement_rule)?;
let sorter = CipSort::from_rules(ref_sort_rules);
let mut queue = vec![beg];
let mut pos = 0_usize;
while pos < queue.len() {
let node = queue[pos];
plist.add(digraph.node(node).get_aux());
let mut edges = digraph.node_edges(node)?;
sorter.prioritize(digraph, context, node, &mut edges, true)?;
for edge in edges {
if digraph.edge(edge).is_beg(node)
&& !digraph.node(digraph.edge(edge).get_end()).is_terminal()
{
queue.push(digraph.edge(edge).get_end());
}
}
pos += 1;
}
Ok(())
}
fn compare_pairs(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipNodeId,
b: CipNodeId,
ref_a: Descriptor,
ref_b: Descriptor,
) -> Result<i32, CipLabelerError> {
let replacement_a = CipRule4b::with_ref(ref_a);
let replacement_b = CipRule4b::with_ref(ref_b);
let a_sorter = CipSort::from_rules(self.get_ref_sorter(sort_rules, &replacement_a)?);
let b_sorter = CipSort::from_rules(self.get_ref_sorter(sort_rules, &replacement_b)?);
let mut a_queue = vec![a];
let mut b_queue = vec![b];
let mut pos = 0_usize;
while pos < a_queue.len() && pos < b_queue.len() {
let a_node = a_queue[pos];
let b_node = b_queue[pos];
let des_a = CipPairList::ref_descriptor(digraph.node(a_node).get_aux());
let des_b = CipPairList::ref_descriptor(digraph.node(b_node).get_aux());
if des_a == ref_a && des_b != ref_b {
return Ok(1);
} else if des_a != ref_a && des_b == ref_b {
return Ok(-1);
}
let mut edges = digraph.node_edges(a_node)?;
a_sorter.prioritize(digraph, context, a_node, &mut edges, true)?;
for edge in edges {
if digraph.edge(edge).is_beg(a_node)
&& !digraph.node(digraph.edge(edge).get_end()).is_terminal()
{
a_queue.push(digraph.edge(edge).get_end());
}
}
let mut edges = digraph.node_edges(b_node)?;
b_sorter.prioritize(digraph, context, b_node, &mut edges, true)?;
for edge in edges {
if digraph.edge(edge).is_beg(b_node)
&& !digraph.node(digraph.edge(edge).get_end()).is_terminal()
{
b_queue.push(digraph.edge(edge).get_end());
}
}
pos += 1;
}
Ok(0)
}
fn get_ref_sorter<'r>(
&'r self,
sort_rules: Option<&'r [&'r dyn CipSequenceRule]>,
replacement_rule: &'r dyn CipSequenceRule,
) -> Result<Vec<&'r dyn CipSequenceRule>, CipLabelerError> {
let mut new_rules = Vec::new();
if let Some(sort_rules) = sort_rules {
new_rules.reserve(sort_rules.len());
let self_ptr = self as &dyn CipSequenceRule as *const dyn CipSequenceRule;
let mut found = false;
for rule in sort_rules {
let rule_ptr = *rule as *const dyn CipSequenceRule;
if std::ptr::addr_eq(rule_ptr, self_ptr) {
found = true;
} else {
new_rules.push(*rule);
}
}
if !found {
return Err(CipLabelerError::Rule4bInstanceNotInRuleSet);
}
}
new_rules.push(replacement_rule);
Ok(new_rules)
}
}
impl CipSequenceRule for CipRule4b {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
self.compare_with_sort_rules(None, digraph, context, a, b)
}
fn compare_with_sort_rules(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_beg = digraph.edge(a).get_beg();
let a_end = digraph.edge(a).get_end();
let b_beg = digraph.edge(b).get_beg();
let b_end = digraph.edge(b).get_end();
if digraph.get_current_root() != a_beg || digraph.get_current_root() != b_beg {
if self.ref_descriptor == Descriptor::None {
return Ok(0);
}
let a_desc = digraph.node(a_end).get_aux();
let b_desc = digraph.node(b_end).get_aux();
if a_desc != Descriptor::None
&& b_desc != Descriptor::None
&& a_desc != Descriptor::ns
&& b_desc != Descriptor::ns
{
let alike = CipPairList::ref_descriptor(self.ref_descriptor)
== CipPairList::ref_descriptor(a_desc);
let blike = CipPairList::ref_descriptor(self.ref_descriptor)
== CipPairList::ref_descriptor(b_desc);
if alike && !blike {
return Ok(1);
}
if blike && !alike {
return Ok(-1);
}
}
return Ok(0);
}
let mut list1 = self
.new_pair_lists(&self.get_reference_descriptors(sort_rules, digraph, context, a_end)?);
let mut list2 = self
.new_pair_lists(&self.get_reference_descriptors(sort_rules, digraph, context, b_end)?);
if list1.is_empty() != list2.is_empty() {
return Err(CipLabelerError::SubstituentsShouldBeTopologicallyEquivalent);
}
if list1.len() == 1 {
self.compare_pairs(
sort_rules,
digraph,
context,
a_end,
b_end,
list1[0].get_ref_descriptor(),
list2[0].get_ref_descriptor(),
)
} else if list1.len() > 1 {
for plist in &mut list1 {
self.fill_pairs(sort_rules, digraph, context, a_end, plist)?;
}
for plist in &mut list2 {
self.fill_pairs(sort_rules, digraph, context, b_end, plist)?;
}
list1.sort_by(|left, right| right.cmp(left));
list2.sort_by(|left, right| right.cmp(left));
for (left, right) in list1.iter().zip(list2.iter()) {
let cmp = left.compare_to(right)?;
if cmp != 0 {
return Ok(cmp);
}
}
Ok(0)
} else {
Ok(0)
}
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule4c;
impl CipRule4c {
fn ord(lab: Descriptor) -> i32 {
match lab {
Descriptor::m | Descriptor::r => 2,
Descriptor::p | Descriptor::s => 1,
_ => 0,
}
}
}
impl CipSequenceRule for CipRule4c {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_ordinal = Self::ord(self.get_bond_label(digraph.edge(a)));
let b_ordinal = Self::ord(self.get_bond_label(digraph.edge(b)));
let cmp = three_way_comparison_i32(a_ordinal, b_ordinal);
if cmp != 0 {
return Ok(cmp);
}
let a_ordinal = Self::ord(digraph.node(digraph.edge(a).get_end()).get_aux());
let b_ordinal = Self::ord(digraph.node(digraph.edge(b).get_end()).get_aux());
Ok(three_way_comparison_i32(a_ordinal, b_ordinal))
}
}
#[derive(Debug, Clone, Copy)]
pub(crate) struct CipRule5New {
ref_descriptor: Descriptor,
}
impl CipRule5New {
pub(crate) fn new() -> Self {
Self {
ref_descriptor: Descriptor::None,
}
}
pub(crate) fn with_ref(ref_descriptor: Descriptor) -> Self {
Self { ref_descriptor }
}
fn fill_pairs(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
beg: CipNodeId,
plist: &mut CipPairList,
) -> Result<(), CipLabelerError> {
let replacement_rule = CipRule5New::with_ref(plist.get_ref_descriptor());
let ref_sort_rules = self.get_ref_sorter(sort_rules, &replacement_rule)?;
let sorter = CipSort::from_rules(ref_sort_rules);
let mut queue = vec![beg];
let mut pos = 0_usize;
while pos < queue.len() {
let node = queue[pos];
plist.add(digraph.node(node).get_aux());
let mut edges = digraph.node_edges(node)?;
sorter.prioritize(digraph, context, node, &mut edges, true)?;
for edge in edges {
if digraph.edge(edge).is_beg(node)
&& !digraph.node(digraph.edge(edge).get_end()).is_terminal()
{
queue.push(digraph.edge(edge).get_end());
}
}
pos += 1;
}
Ok(())
}
fn get_ref_sorter<'r>(
&'r self,
sort_rules: Option<&'r [&'r dyn CipSequenceRule]>,
replacement_rule: &'r dyn CipSequenceRule,
) -> Result<Vec<&'r dyn CipSequenceRule>, CipLabelerError> {
let mut new_rules = Vec::new();
if let Some(sort_rules) = sort_rules {
new_rules.reserve(sort_rules.len());
let self_ptr = self as &dyn CipSequenceRule as *const dyn CipSequenceRule;
let mut found = false;
for rule in sort_rules {
let rule_ptr = *rule as *const dyn CipSequenceRule;
if std::ptr::addr_eq(rule_ptr, self_ptr) {
found = true;
} else {
new_rules.push(*rule);
}
}
if !found {
return Err(CipLabelerError::Rule5NewInstanceNotInRuleSet);
}
}
new_rules.push(replacement_rule);
Ok(new_rules)
}
}
impl CipSequenceRule for CipRule5New {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
self.compare_with_sort_rules(None, digraph, context, a, b)
}
fn compare_with_sort_rules(
&self,
sort_rules: Option<&[&dyn CipSequenceRule]>,
digraph: &mut CipDigraph<'_>,
context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_beg = digraph.edge(a).get_beg();
let a_end = digraph.edge(a).get_end();
let b_beg = digraph.edge(b).get_beg();
let b_end = digraph.edge(b).get_end();
if digraph.get_current_root() != a_beg || digraph.get_current_root() != b_beg {
if self.ref_descriptor == Descriptor::None {
return Ok(0);
}
let a_desc = digraph.node(a_end).get_aux();
let b_desc = digraph.node(b_end).get_aux();
if a_desc != Descriptor::None
&& b_desc != Descriptor::None
&& a_desc != Descriptor::ns
&& b_desc != Descriptor::ns
{
let alike = CipPairList::ref_descriptor(self.ref_descriptor)
== CipPairList::ref_descriptor(a_desc);
let blike = CipPairList::ref_descriptor(self.ref_descriptor)
== CipPairList::ref_descriptor(b_desc);
if alike && !blike {
return Ok(1);
}
if blike && !alike {
return Ok(-1);
}
}
return Ok(0);
}
let mut list_ra = CipPairList::with_ref(Descriptor::R);
let mut list_rb = CipPairList::with_ref(Descriptor::R);
let mut list_sa = CipPairList::with_ref(Descriptor::S);
let mut list_sb = CipPairList::with_ref(Descriptor::S);
self.fill_pairs(sort_rules, digraph, context, a_end, &mut list_ra)?;
self.fill_pairs(sort_rules, digraph, context, a_end, &mut list_sa)?;
self.fill_pairs(sort_rules, digraph, context, b_end, &mut list_rb)?;
self.fill_pairs(sort_rules, digraph, context, b_end, &mut list_sb)?;
let cmp_r = list_ra.compare_to(&list_rb)?;
let cmp_s = list_sa.compare_to(&list_sb)?;
if cmp_r < 0 {
Ok(if cmp_s < 0 { -1 } else { -2 })
} else if cmp_r > 0 {
Ok(if cmp_s > 0 { 1 } else { 2 })
} else {
Ok(0)
}
}
}
#[derive(Debug, Default, Clone, Copy)]
pub(crate) struct CipRule6;
impl CipSequenceRule for CipRule6 {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let Some(ref_atom) = digraph.get_rule6_ref() else {
return Ok(0);
};
let a_atom = digraph.node(digraph.edge(a).get_end()).atom_idx();
let b_atom = digraph.node(digraph.edge(b).get_end()).atom_idx();
if Some(ref_atom) == a_atom && Some(ref_atom) != b_atom {
Ok(1)
} else if Some(ref_atom) != a_atom && Some(ref_atom) == b_atom {
Ok(-1)
} else {
Ok(0)
}
}
}
pub(crate) struct CipConfiguration<'a> {
foci: Vec<usize>,
carriers: Vec<Option<usize>>,
digraph: CipDigraph<'a>,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) struct CipAtomPrimaryLabel {
atom_idx: usize,
cip_code: &'static str,
cip_neighbor_order: Vec<usize>,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) struct CipBondPrimaryLabel {
bond_idx: usize,
stereo_atoms: [usize; 2],
stereo: BondStereo,
cip_code: &'static str,
cip_neighbor_order: Vec<usize>,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) struct CipAtropisomerBondPrimaryLabel {
bond_idx: usize,
cip_code: &'static str,
cip_neighbor_order: Vec<usize>,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) enum CipPrimaryLabel {
Atom(CipAtomPrimaryLabel),
Bond(CipBondPrimaryLabel),
AtropisomerBond(CipAtropisomerBondPrimaryLabel),
}
pub(crate) enum CipConfig<'a> {
Tetrahedral(CipTetrahedral<'a>),
Sp2Bond(CipSp2Bond<'a>),
AtropisomerBond(CipAtropisomerBond<'a>),
}
impl<'a> CipConfig<'a> {
fn get_foci(&self) -> &[usize] {
match self {
Self::Tetrahedral(config) => config.get_foci(),
Self::Sp2Bond(config) => config.get_foci(),
Self::AtropisomerBond(config) => config.get_foci(),
}
}
fn get_digraph_mut(&mut self) -> &mut CipDigraph<'a> {
match self {
Self::Tetrahedral(config) => config.configuration.get_digraph(),
Self::Sp2Bond(config) => config.configuration.get_digraph(),
Self::AtropisomerBond(config) => config.configuration.get_digraph(),
}
}
fn reset_primary_label(&mut self) {
match self {
Self::Tetrahedral(config) => config.reset_primary_label(),
Self::Sp2Bond(config) => config.reset_primary_label(),
Self::AtropisomerBond(config) => config.reset_primary_label(),
}
}
fn has_primary_label(&self) -> bool {
match self {
Self::Tetrahedral(config) => config.has_primary_label(),
Self::Sp2Bond(config) => config.has_primary_label(),
Self::AtropisomerBond(config) => config.has_primary_label(),
}
}
fn label(
&mut self,
rules: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
match self {
Self::Tetrahedral(config) => config.label(rules, context),
Self::Sp2Bond(config) => config.label(rules, context),
Self::AtropisomerBond(config) => config.label(rules, context),
}
}
fn label_with_external_digraph(
&mut self,
node: CipNodeId,
digraph: &mut CipDigraph<'_>,
rules: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
match self {
Self::Tetrahedral(config) => {
config.label_with_external_digraph(node, digraph, rules, context)
}
Self::Sp2Bond(config) => {
config.label_with_external_digraph(node, digraph, rules, context)
}
Self::AtropisomerBond(config) => {
config.label_with_external_digraph(node, digraph, rules, context)
}
}
}
fn set_primary_label(&mut self, desc: Descriptor) -> Result<(), CipLabelerError> {
match self {
Self::Tetrahedral(config) => config.set_primary_label(desc),
Self::Sp2Bond(config) => config.set_primary_label(desc),
Self::AtropisomerBond(config) => config.set_primary_label(desc),
}
}
fn primary_label(&self) -> Option<CipPrimaryLabel> {
match self {
Self::Tetrahedral(config) => config.primary_label().cloned().map(CipPrimaryLabel::Atom),
Self::Sp2Bond(config) => config.primary_label().cloned().map(CipPrimaryLabel::Bond),
Self::AtropisomerBond(config) => config
.primary_label()
.cloned()
.map(CipPrimaryLabel::AtropisomerBond),
}
}
}
impl<'a> CipConfiguration<'a> {
pub(crate) const IMPLICIT_H: usize = usize::MAX;
pub(crate) fn parity4<T: PartialEq>(
trg: &[T],
reference: &[T],
) -> Result<i32, CipLabelerError> {
if reference.len() != 4 || trg.len() != reference.len() {
return Err(CipLabelerError::ParityVectorsMustHaveSize4);
}
let r = reference;
let t = trg;
if r[0] == t[0] {
if r[1] == t[1] {
if r[2] == t[2] && r[3] == t[3] {
return Ok(2);
}
if r[2] == t[3] && r[3] == t[2] {
return Ok(1);
}
} else if r[1] == t[2] {
if r[2] == t[1] && r[3] == t[3] {
return Ok(1);
}
if r[2] == t[3] && r[3] == t[1] {
return Ok(2);
}
} else if r[1] == t[3] {
if r[2] == t[2] && r[3] == t[1] {
return Ok(1);
}
if r[2] == t[1] && r[3] == t[2] {
return Ok(2);
}
}
} else if r[0] == t[1] {
if r[1] == t[0] {
if r[2] == t[2] && r[3] == t[3] {
return Ok(1);
}
if r[2] == t[3] && r[3] == t[2] {
return Ok(2);
}
} else if r[1] == t[2] {
if r[2] == t[0] && r[3] == t[3] {
return Ok(2);
}
if r[2] == t[3] && r[3] == t[0] {
return Ok(1);
}
} else if r[1] == t[3] {
if r[2] == t[2] && r[3] == t[0] {
return Ok(2);
}
if r[2] == t[0] && r[3] == t[2] {
return Ok(1);
}
}
} else if r[0] == t[2] {
if r[1] == t[1] {
if r[2] == t[0] && r[3] == t[3] {
return Ok(1);
}
if r[2] == t[3] && r[3] == t[0] {
return Ok(2);
}
} else if r[1] == t[0] {
if r[2] == t[1] && r[3] == t[3] {
return Ok(2);
}
if r[2] == t[3] && r[3] == t[1] {
return Ok(1);
}
} else if r[1] == t[3] {
if r[2] == t[0] && r[3] == t[1] {
return Ok(2);
}
if r[2] == t[1] && r[3] == t[0] {
return Ok(1);
}
}
} else if r[0] == t[3] {
if r[1] == t[1] {
if r[2] == t[2] && r[3] == t[0] {
return Ok(1);
}
if r[2] == t[0] && r[3] == t[2] {
return Ok(2);
}
} else if r[1] == t[2] {
if r[2] == t[1] && r[3] == t[0] {
return Ok(2);
}
if r[2] == t[0] && r[3] == t[1] {
return Ok(1);
}
} else if r[1] == t[0] {
if r[2] == t[2] && r[3] == t[1] {
return Ok(2);
}
if r[2] == t[1] && r[3] == t[2] {
return Ok(1);
}
}
}
Ok(0)
}
pub(crate) fn new(molecule: &'a Molecule, focus: usize) -> Result<Self, CipLabelerError> {
Ok(Self {
foci: vec![focus],
carriers: Vec::new(),
digraph: CipDigraph::new(molecule, focus, false)?,
})
}
pub(crate) fn with_foci(
molecule: &'a Molecule,
foci: Vec<usize>,
atropisomer_mode: bool,
) -> Result<Self, CipLabelerError> {
let focus = *foci
.first()
.ok_or(CipLabelerError::EmptyConfigurationFoci)?;
Ok(Self {
foci,
carriers: Vec::new(),
digraph: CipDigraph::new(molecule, focus, atropisomer_mode)?,
})
}
pub(crate) fn get_focus(&self) -> usize {
self.foci[0]
}
pub(crate) fn get_foci(&self) -> &[usize] {
&self.foci
}
pub(crate) fn set_carriers(&mut self, carriers: Vec<Option<usize>>) {
self.carriers = carriers;
}
pub(crate) fn get_carriers(&self) -> &[Option<usize>] {
&self.carriers
}
pub(crate) fn get_digraph(&mut self) -> &mut CipDigraph<'a> {
&mut self.digraph
}
pub(crate) fn label(
&self,
_node: CipNodeId,
_digraph: &mut CipDigraph<'_>,
_comp: &CipRules,
) -> Descriptor {
Descriptor::Unknown
}
fn find_internal_edge(
digraph: &CipDigraph<'_>,
edges: &[CipEdgeId],
f1: usize,
f2: usize,
) -> Option<CipEdgeId> {
edges.iter().copied().find(|edge_id| {
let edge = digraph.edge(*edge_id);
if digraph.node(edge.get_beg()).is_duplicate()
|| digraph.node(edge.get_end()).is_duplicate()
{
return false;
}
Self::is_internal_edge(digraph, *edge_id, f1, f2)
})
}
fn is_internal_edge(
digraph: &CipDigraph<'_>,
edge_id: CipEdgeId,
f1: usize,
f2: usize,
) -> bool {
let edge = digraph.edge(edge_id);
let beg = digraph.node(edge.get_beg()).atom_idx();
let end = digraph.node(edge.get_end()).atom_idx();
(beg == Some(f1) && end == Some(f2)) || (beg == Some(f2) && end == Some(f1))
}
fn remove_internal_edges(
digraph: &CipDigraph<'_>,
edges: &mut Vec<CipEdgeId>,
f1: usize,
f2: usize,
) {
edges.retain(|edge_id| !Self::is_internal_edge(digraph, *edge_id, f1, f2));
}
fn is_duplicate_or_hydrogen_edge(digraph: &CipDigraph<'_>, edge_id: CipEdgeId) -> bool {
let edge = digraph.edge(edge_id);
digraph.node(edge.get_beg()).is_duplicate_or_h()
|| digraph.node(edge.get_end()).is_duplicate_or_h()
}
fn remove_duplicates_and_hs(digraph: &CipDigraph<'_>, edges: &mut Vec<CipEdgeId>) {
edges.retain(|edge_id| !Self::is_duplicate_or_hydrogen_edge(digraph, *edge_id));
}
}
pub(crate) struct CipTetrahedral<'a> {
configuration: CipConfiguration<'a>,
ranked_anchors: Vec<usize>,
primary_label: Option<CipAtomPrimaryLabel>,
}
impl<'a> CipTetrahedral<'a> {
pub(crate) fn new(molecule: &'a Molecule, focus: usize) -> Result<Self, CipLabelerError> {
let mut configuration = CipConfiguration::new(molecule, focus)?;
let atom = configuration.digraph.mol().atom(focus)?;
match atom.chiral_tag() {
ChiralTag::TetrahedralCcw | ChiralTag::TetrahedralCw => {}
ChiralTag::Unspecified
| ChiralTag::Other
| ChiralTag::Tetrahedral
| ChiralTag::Allene
| ChiralTag::SquarePlanar
| ChiralTag::TrigonalBipyramidal
| ChiralTag::Octahedral => return Err(CipLabelerError::BadTetrahedralConfig),
}
let mut carriers = Vec::with_capacity(4);
for nbr in configuration.digraph.mol().neighbor_indices(focus)? {
carriers.push(Some(nbr));
}
if carriers.len() < 4 {
carriers.push(Some(focus));
}
if carriers.len() < 4 {
carriers.push(None);
}
if carriers.len() != 4 {
return Err(CipLabelerError::TetrahedralConfigurationMustHave4Carriers);
}
configuration.set_carriers(carriers);
Ok(Self {
configuration,
ranked_anchors: Vec::new(),
primary_label: None,
})
}
pub(crate) fn get_focus(&self) -> usize {
self.configuration.get_focus()
}
pub(crate) fn get_foci(&self) -> &[usize] {
self.configuration.get_foci()
}
pub(crate) fn get_carriers(&self) -> &[Option<usize>] {
self.configuration.get_carriers()
}
pub(crate) fn ranked_anchors(&self) -> &[usize] {
&self.ranked_anchors
}
pub(crate) fn primary_label(&self) -> Option<&CipAtomPrimaryLabel> {
self.primary_label.as_ref()
}
pub(crate) fn set_primary_label(&mut self, desc: Descriptor) -> Result<(), CipLabelerError> {
match desc {
Descriptor::R | Descriptor::S | Descriptor::r | Descriptor::s => {
self.primary_label = Some(CipAtomPrimaryLabel {
atom_idx: self.configuration.get_focus(),
cip_code: descriptor_to_string(desc),
cip_neighbor_order: self.ranked_anchors.clone(),
});
Ok(())
}
Descriptor::seqTrans
| Descriptor::seqCis
| Descriptor::E
| Descriptor::Z
| Descriptor::M
| Descriptor::P
| Descriptor::m
| Descriptor::p
| Descriptor::SP_4
| Descriptor::TBPY_5
| Descriptor::OC_6 => Err(CipLabelerError::DescriptorNotSupportedForAtoms),
Descriptor::None | Descriptor::Unknown | Descriptor::ns => {
Err(CipLabelerError::InvalidAtomDescriptor)
}
}
}
pub(crate) fn has_primary_label(&self) -> bool {
self.primary_label.is_some()
|| self
.configuration
.digraph
.mol()
.atom(self.configuration.get_focus())
.is_ok_and(|atom| atom.prop("_CIPCode").is_some())
}
pub(crate) fn reset_primary_label(&mut self) {
self.primary_label = None;
}
pub(crate) fn label(
&mut self,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let root = self.configuration.digraph.get_original_root();
if self.configuration.digraph.get_current_root() != root {
self.configuration.digraph.change_root(root)?;
}
self.label_node(root, comp, context)
}
pub(crate) fn label_with_external_digraph(
&mut self,
node: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
digraph.change_root(node)?;
self.label_node_in_digraph(node, digraph, comp, context)
}
fn label_node(
&mut self,
node: CipNodeId,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let focus = self.configuration.get_focus();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
focus,
&carriers,
&mut self.ranked_anchors,
node,
&mut self.configuration.digraph,
comp,
context,
)
}
fn label_node_in_digraph(
&mut self,
node: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let focus = self.configuration.get_focus();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
focus,
&carriers,
&mut self.ranked_anchors,
node,
digraph,
comp,
context,
)
}
fn label_node_impl(
focus: usize,
carriers: &[Option<usize>],
ranked_anchors: &mut Vec<usize>,
node: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let mut edges = digraph.node_edges(node)?;
ranked_anchors.clear();
if edges.len() < 3 {
return Ok(Descriptor::ns);
}
let mut priority = comp.sort(digraph, context, node, &mut edges, true)?;
let is_unique = priority.is_unique();
if !is_unique && edges.len() == 4 {
if comp.get_num_sub_rules() == 3 {
return Ok(Descriptor::Unknown);
}
let partition = comp.get_sorter().get_groups(digraph, context, &edges)?;
if partition.len() == 2 {
let ref_atom = digraph.node(digraph.edge(edges[1]).get_end()).atom_idx();
digraph.set_rule6_ref(ref_atom)?;
priority = comp.sort(digraph, context, node, &mut edges, true)?;
digraph.set_rule6_ref(None)?;
} else if partition.len() == 1 {
let ref_atom = digraph.node(digraph.edge(edges[0]).get_end()).atom_idx();
digraph.set_rule6_ref(ref_atom)?;
comp.sort(digraph, context, node, &mut edges, true)?;
let nbrs1 = edges.clone();
let ref_atom = digraph.node(digraph.edge(edges[1]).get_end()).atom_idx();
digraph.set_rule6_ref(ref_atom)?;
priority = comp.sort(digraph, context, node, &mut edges, true)?;
digraph.set_rule6_ref(None)?;
if CipConfiguration::parity4(&nbrs1, &edges)? == 1 {
return Ok(Descriptor::Unknown);
}
}
if !priority.is_unique() {
return Ok(Descriptor::Unknown);
}
} else if !is_unique {
return Ok(Descriptor::Unknown);
}
let mut ordered = vec![None; 4];
let mut idx = 0_usize;
ranked_anchors.reserve(4);
for edge in &edges {
let end = digraph.edge(*edge).get_end();
if digraph.node(end).is_set(CipNode::BOND_DUPLICATE)
|| digraph.node(end).is_set(CipNode::IMPL_HYDROGEN)
{
continue;
}
let atom = digraph.node(end).atom_idx();
if idx < 4 {
ordered[idx] = atom;
}
if let Some(atom_idx) = atom {
ranked_anchors.push(atom_idx);
}
idx += 1;
}
if idx < 4 {
ordered[idx] = Some(focus);
}
let parity = CipConfiguration::parity4(&ordered, carriers)?;
if parity == 0 {
return Err(CipLabelerError::CarrierMismatch);
}
let mut config = digraph.mol().atom(focus)?.chiral_tag();
if parity == 1 {
config = match config {
ChiralTag::TetrahedralCcw => ChiralTag::TetrahedralCw,
ChiralTag::TetrahedralCw => ChiralTag::TetrahedralCcw,
_ => config,
};
}
if config == ChiralTag::TetrahedralCcw {
if priority.is_pseudo_asymetric() {
Ok(Descriptor::s)
} else {
Ok(Descriptor::S)
}
} else if config == ChiralTag::TetrahedralCw {
if priority.is_pseudo_asymetric() {
Ok(Descriptor::r)
} else {
Ok(Descriptor::R)
}
} else {
Ok(Descriptor::Unknown)
}
}
}
pub(crate) struct CipSp2Bond<'a> {
configuration: CipConfiguration<'a>,
bond_idx: usize,
cfg: BondStereo,
ranked_anchors: Vec<usize>,
primary_label: Option<CipBondPrimaryLabel>,
}
impl<'a> CipSp2Bond<'a> {
pub(crate) fn new(
molecule: &'a Molecule,
bond_idx: usize,
start_atom: usize,
end_atom: usize,
cfg: BondStereo,
) -> Result<Self, CipLabelerError> {
let mut configuration =
CipConfiguration::with_foci(molecule, vec![start_atom, end_atom], false)?;
configuration.digraph.mol().atom(start_atom)?;
configuration.digraph.mol().atom(end_atom)?;
let bond = configuration.digraph.mol().bond(bond_idx)?;
if bond.order() != BondOrder::Double {
return Err(CipLabelerError::BadSp2BondFoci);
}
if !((bond.begin().index() == start_atom && bond.end().index() == end_atom)
|| (bond.begin().index() == end_atom && bond.end().index() == start_atom))
{
return Err(CipLabelerError::BadSp2BondFoci);
}
if !matches!(cfg, BondStereo::Trans | BondStereo::Cis) {
return Err(CipLabelerError::BadSp2BondConfig);
}
let stereo_atoms = if let Some([left, right]) = bond.stereo_atoms() {
[left.index(), right.index()]
} else if matches!(bond.stereo(), BondStereo::E | BondStereo::Z) {
let left = Self::find_highest_cip_neighbor_like_rdkit(
configuration.digraph.mol(),
start_atom,
end_atom,
)?;
let right = Self::find_highest_cip_neighbor_like_rdkit(
configuration.digraph.mol(),
end_atom,
start_atom,
)?;
match (left, right) {
(Some(left), Some(right)) => [left, right],
_ => return Err(CipLabelerError::IncorrectNumberOfStereoAtoms),
}
} else {
return Err(CipLabelerError::IncorrectNumberOfStereoAtoms);
};
configuration.digraph.mol().atom(stereo_atoms[0])?;
configuration.digraph.mol().atom(stereo_atoms[1])?;
configuration.set_carriers(vec![Some(stereo_atoms[0]), Some(stereo_atoms[1])]);
Ok(Self {
configuration,
bond_idx,
cfg,
ranked_anchors: Vec::new(),
primary_label: None,
})
}
fn find_highest_cip_neighbor_like_rdkit(
mol: &CipMol<'_>,
atom_idx: usize,
skip_atom_idx: usize,
) -> Result<Option<usize>, CipLabelerError> {
mol.atom(atom_idx)?;
let mut best_cip_rank = 0_u32;
let mut best_cip_ranked_atom = None;
for neighbor in mol.neighbor_indices(atom_idx)? {
if neighbor == skip_atom_idx {
continue;
}
let Some(rank_text) = mol.atom(neighbor)?.prop("_CIPRank") else {
return Ok(None);
};
let Ok(cip) = rank_text.parse::<u32>() else {
return Ok(None);
};
if cip > best_cip_rank || best_cip_ranked_atom.is_none() {
best_cip_rank = cip;
best_cip_ranked_atom = Some(neighbor);
} else if cip == best_cip_rank {
best_cip_ranked_atom = None;
}
}
Ok(best_cip_ranked_atom)
}
pub(crate) fn get_foci(&self) -> &[usize] {
self.configuration.get_foci()
}
pub(crate) fn get_carriers(&self) -> &[Option<usize>] {
self.configuration.get_carriers()
}
pub(crate) fn ranked_anchors(&self) -> &[usize] {
&self.ranked_anchors
}
pub(crate) fn primary_label(&self) -> Option<&CipBondPrimaryLabel> {
self.primary_label.as_ref()
}
pub(crate) fn set_primary_label(&mut self, desc: Descriptor) -> Result<(), CipLabelerError> {
match desc {
Descriptor::seqTrans | Descriptor::E | Descriptor::seqCis | Descriptor::Z => {
let carriers = self.configuration.get_carriers();
let stereo_atoms = [
carriers
.first()
.and_then(|carrier| *carrier)
.ok_or(CipLabelerError::IncorrectNumberOfStereoAtoms)?,
carriers
.get(1)
.and_then(|carrier| *carrier)
.ok_or(CipLabelerError::IncorrectNumberOfStereoAtoms)?,
];
self.primary_label = Some(CipBondPrimaryLabel {
bond_idx: self.bond_idx,
stereo_atoms,
stereo: self.cfg,
cip_code: descriptor_to_string(desc),
cip_neighbor_order: self.ranked_anchors.clone(),
});
Ok(())
}
Descriptor::R
| Descriptor::S
| Descriptor::r
| Descriptor::s
| Descriptor::M
| Descriptor::P
| Descriptor::m
| Descriptor::p
| Descriptor::SP_4
| Descriptor::TBPY_5
| Descriptor::OC_6 => Err(CipLabelerError::DescriptorNotSupportedForDoubleBonds),
Descriptor::None | Descriptor::Unknown | Descriptor::ns => {
Err(CipLabelerError::InvalidBondDescriptor)
}
}
}
pub(crate) fn has_primary_label(&self) -> bool {
self.primary_label.is_some()
|| self
.configuration
.digraph
.mol()
.bond(self.bond_idx)
.is_ok_and(|bond| bond.prop("_CIPCode").is_some())
}
pub(crate) fn reset_primary_label(&mut self) {
self.primary_label = None;
}
pub(crate) fn label(
&mut self,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let root1 = self.configuration.digraph.get_original_root();
if self.configuration.digraph.get_current_root() != root1 {
self.configuration.digraph.change_root(root1)?;
}
self.label_node(root1, comp, context)
}
pub(crate) fn label_with_external_digraph(
&mut self,
root1: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let foci = self.configuration.get_foci().to_vec();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
&foci,
&carriers,
self.cfg,
&mut self.ranked_anchors,
root1,
digraph,
comp,
context,
)
}
fn label_node(
&mut self,
root1: CipNodeId,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let foci = self.configuration.get_foci().to_vec();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
&foci,
&carriers,
self.cfg,
&mut self.ranked_anchors,
root1,
&mut self.configuration.digraph,
comp,
context,
)
}
fn label_node_impl(
foci: &[usize],
carriers: &[Option<usize>],
cfg: BondStereo,
ranked_anchors: &mut Vec<usize>,
root1: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let focus1 = foci[0];
let focus2 = foci[1];
ranked_anchors.clear();
let root1_edges = digraph.node_edges(root1)?;
let Some(internal) =
CipConfiguration::find_internal_edge(digraph, &root1_edges, focus1, focus2)
else {
return Ok(Descriptor::Unknown);
};
let root2 = digraph.edge(internal).get_other(internal, root1)?;
let mut edges1 = digraph.node_edges(root1)?;
let mut edges2 = digraph.node_edges(root2)?;
CipConfiguration::remove_internal_edges(digraph, &mut edges1, focus1, focus2);
CipConfiguration::remove_internal_edges(digraph, &mut edges2, focus1, focus2);
if edges1.is_empty() || edges2.is_empty() {
return Ok(Descriptor::Unknown);
}
let mut carriers = carriers.to_vec();
let mut config = cfg;
if digraph.node(root1).atom_idx() == Some(focus2) {
carriers.swap(0, 1);
}
digraph.change_root(root1)?;
let priority1 = comp.sort(digraph, context, root1, &mut edges1, true)?;
if !priority1.is_unique() {
return Ok(Descriptor::Unknown);
}
if edges1.len() > 1
&& carriers[0] != digraph.node(digraph.edge(edges1[0]).get_end()).atom_idx()
{
config = match config {
BondStereo::Cis => BondStereo::Trans,
BondStereo::Trans => BondStereo::Cis,
_ => config,
};
}
digraph.change_root(root2)?;
let priority2 = comp.sort(digraph, context, root2, &mut edges2, true)?;
if !priority2.is_unique() {
return Ok(Descriptor::Unknown);
}
if edges2.len() > 1
&& carriers[1] != digraph.node(digraph.edge(edges2[0]).get_end()).atom_idx()
{
config = match config {
BondStereo::Cis => BondStereo::Trans,
BondStereo::Trans => BondStereo::Cis,
_ => config,
};
}
let carrier1_idx = digraph
.node(digraph.edge(edges1[0]).get_end())
.atom_idx()
.unwrap_or(CipNode::NO_ATOM_INDEX);
let carrier2_idx = digraph
.node(digraph.edge(edges2[0]).get_end())
.atom_idx()
.unwrap_or(CipNode::NO_ATOM_INDEX);
if digraph.node(digraph.edge(edges1[0]).get_beg()).atom_idx() == Some(focus1) {
ranked_anchors.extend([carrier1_idx, carrier2_idx]);
} else if digraph.node(digraph.edge(edges2[0]).get_beg()).atom_idx() == Some(focus1) {
ranked_anchors.extend([carrier2_idx, carrier1_idx]);
}
if config == BondStereo::Cis {
if priority1.is_pseudo_asymetric() != priority2.is_pseudo_asymetric() {
Ok(Descriptor::seqCis)
} else {
Ok(Descriptor::Z)
}
} else if config == BondStereo::Trans {
if priority1.is_pseudo_asymetric() != priority2.is_pseudo_asymetric() {
Ok(Descriptor::seqTrans)
} else {
Ok(Descriptor::E)
}
} else {
Ok(Descriptor::Unknown)
}
}
}
pub(crate) struct CipAtropisomerBond<'a> {
configuration: CipConfiguration<'a>,
bond_idx: usize,
cfg: BondStereo,
ranked_anchors: Vec<usize>,
primary_label: Option<CipAtropisomerBondPrimaryLabel>,
}
impl<'a> CipAtropisomerBond<'a> {
pub(crate) fn new(
molecule: &'a Molecule,
bond_idx: usize,
start_atom: usize,
end_atom: usize,
cfg: BondStereo,
) -> Result<Self, CipLabelerError> {
let mut configuration =
CipConfiguration::with_foci(molecule, vec![start_atom, end_atom], true)?;
configuration.digraph.mol().atom(start_atom)?;
configuration.digraph.mol().atom(end_atom)?;
let bond = configuration.digraph.mol().bond(bond_idx)?;
if !((bond.begin().index() == start_atom && bond.end().index() == end_atom)
|| (bond.begin().index() == end_atom && bond.end().index() == start_atom))
{
return Err(CipLabelerError::BadAtropisomerBondFoci);
}
if !matches!(cfg, BondStereo::AtropCw | BondStereo::AtropCcw) {
return Err(CipLabelerError::BadAtropisomerBondConfig);
}
if let Some(carriers) =
Self::atropisomer_carriers_like_rdkit(configuration.digraph.mol(), bond_idx)?
{
configuration.set_carriers(vec![Some(carriers[0]), Some(carriers[1])]);
}
Ok(Self {
configuration,
bond_idx,
cfg,
ranked_anchors: Vec::new(),
primary_label: None,
})
}
fn atropisomer_carriers_like_rdkit(
mol: &CipMol<'_>,
bond_idx: usize,
) -> Result<Option<[usize; 2]>, CipLabelerError> {
let bond = mol.bond(bond_idx)?;
let foci = [bond.begin().index(), bond.end().index()];
let mut carriers = [0_usize; 2];
for (side, focus) in foci.into_iter().enumerate() {
let mut nbr_bonds = mol
.bond_indices_for_atom(focus)?
.into_iter()
.filter(|candidate| *candidate != bond_idx)
.collect::<Vec<_>>();
if nbr_bonds.is_empty() {
return Ok(None);
}
if nbr_bonds.len() == 2 {
let other0 = mol.other_atom_idx(nbr_bonds[0], focus)?;
let other1 = mol.other_atom_idx(nbr_bonds[1], focus)?;
if other1 < other0 {
nbr_bonds.swap(0, 1);
}
}
carriers[side] = mol.other_atom_idx(nbr_bonds[0], focus)?;
}
Ok(Some(carriers))
}
pub(crate) fn get_foci(&self) -> &[usize] {
self.configuration.get_foci()
}
pub(crate) fn get_carriers(&self) -> &[Option<usize>] {
self.configuration.get_carriers()
}
pub(crate) fn ranked_anchors(&self) -> &[usize] {
&self.ranked_anchors
}
pub(crate) fn primary_label(&self) -> Option<&CipAtropisomerBondPrimaryLabel> {
self.primary_label.as_ref()
}
pub(crate) fn set_primary_label(&mut self, desc: Descriptor) -> Result<(), CipLabelerError> {
match desc {
Descriptor::M | Descriptor::P | Descriptor::m | Descriptor::p => {
self.primary_label = Some(CipAtropisomerBondPrimaryLabel {
bond_idx: self.bond_idx,
cip_code: descriptor_to_string(desc),
cip_neighbor_order: self.ranked_anchors.clone(),
});
Ok(())
}
Descriptor::R
| Descriptor::S
| Descriptor::r
| Descriptor::s
| Descriptor::SP_4
| Descriptor::TBPY_5
| Descriptor::OC_6
| Descriptor::seqTrans
| Descriptor::E
| Descriptor::seqCis
| Descriptor::Z => Err(CipLabelerError::DescriptorNotSupportedForAtropisomerBonds),
Descriptor::None | Descriptor::Unknown | Descriptor::ns => {
Err(CipLabelerError::InvalidBondDescriptor)
}
}
}
pub(crate) fn has_primary_label(&self) -> bool {
self.primary_label.is_some()
|| self
.configuration
.digraph
.mol()
.bond(self.bond_idx)
.is_ok_and(|bond| bond.prop("_CIPCode").is_some())
}
pub(crate) fn reset_primary_label(&mut self) {
self.primary_label = None;
}
pub(crate) fn label(
&mut self,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let root1 = self.configuration.digraph.get_original_root();
if self.configuration.digraph.get_current_root() != root1 {
self.configuration.digraph.change_root(root1)?;
}
self.label_with_root(root1, comp, context)
}
pub(crate) fn label_with_external_digraph(
&mut self,
root1: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let foci = self.configuration.get_foci().to_vec();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
&foci,
&carriers,
self.cfg,
&mut self.ranked_anchors,
root1,
digraph,
comp,
context,
)
}
fn label_with_root(
&mut self,
root1: CipNodeId,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let foci = self.configuration.get_foci().to_vec();
let carriers = self.configuration.get_carriers().to_vec();
Self::label_node_impl(
&foci,
&carriers,
self.cfg,
&mut self.ranked_anchors,
root1,
&mut self.configuration.digraph,
comp,
context,
)
}
fn label_node_impl(
foci: &[usize],
carriers: &[Option<usize>],
cfg: BondStereo,
ranked_anchors: &mut Vec<usize>,
root1: CipNodeId,
digraph: &mut CipDigraph<'_>,
comp: &CipRules,
context: &mut CipLabelerContext,
) -> Result<Descriptor, CipLabelerError> {
let focus1 = foci[0];
let focus2 = foci[1];
ranked_anchors.clear();
if carriers.len() < 2 {
return Ok(Descriptor::Unknown);
}
let root1_edges = digraph.node_edges(root1)?;
let Some(internal) =
CipConfiguration::find_internal_edge(digraph, &root1_edges, focus1, focus2)
else {
return Ok(Descriptor::Unknown);
};
let root2 = digraph.edge(internal).get_other(internal, root1)?;
let mut edges1 = digraph.node_edges(root1)?;
let mut edges2 = digraph.node_edges(root2)?;
CipConfiguration::remove_internal_edges(digraph, &mut edges1, focus1, focus2);
CipConfiguration::remove_internal_edges(digraph, &mut edges2, focus1, focus2);
CipConfiguration::remove_duplicates_and_hs(digraph, &mut edges1);
CipConfiguration::remove_duplicates_and_hs(digraph, &mut edges2);
if edges1.is_empty() || edges2.is_empty() {
return Ok(Descriptor::Unknown);
}
let mut carriers = carriers.to_vec();
let mut config = cfg;
if digraph.node(root1).atom_idx() == Some(focus2) {
carriers.swap(0, 1);
}
digraph.change_root(root1)?;
let priority1 = comp.sort(digraph, context, root1, &mut edges1, true)?;
if !priority1.is_unique() {
return Ok(Descriptor::Unknown);
}
if edges1.len() > 1
&& carriers[0] == digraph.node(digraph.edge(edges1[1]).get_end()).atom_idx()
{
config = match config {
BondStereo::AtropCcw => BondStereo::AtropCw,
BondStereo::AtropCw => BondStereo::AtropCcw,
_ => config,
};
}
digraph.change_root(root2)?;
let priority2 = comp.sort(digraph, context, root2, &mut edges2, true)?;
if !priority2.is_unique() {
return Ok(Descriptor::Unknown);
}
if edges2.len() > 1
&& carriers[1] == digraph.node(digraph.edge(edges2[1]).get_end()).atom_idx()
{
config = match config {
BondStereo::AtropCcw => BondStereo::AtropCw,
BondStereo::AtropCw => BondStereo::AtropCcw,
_ => config,
};
}
let carrier1_idx = digraph
.node(digraph.edge(edges1[0]).get_end())
.atom_idx()
.unwrap_or(CipNode::NO_ATOM_INDEX);
let carrier2_idx = digraph
.node(digraph.edge(edges2[0]).get_end())
.atom_idx()
.unwrap_or(CipNode::NO_ATOM_INDEX);
if digraph.node(digraph.edge(edges1[0]).get_beg()).atom_idx() == Some(focus1) {
ranked_anchors.extend([carrier1_idx, carrier2_idx]);
} else if digraph.node(digraph.edge(edges2[0]).get_beg()).atom_idx() == Some(focus1) {
ranked_anchors.extend([carrier2_idx, carrier1_idx]);
}
if config == BondStereo::AtropCcw {
if priority1.is_pseudo_asymetric() || priority2.is_pseudo_asymetric() {
Ok(Descriptor::m)
} else {
Ok(Descriptor::M)
}
} else if config == BondStereo::AtropCw {
if priority1.is_pseudo_asymetric() || priority2.is_pseudo_asymetric() {
Ok(Descriptor::p)
} else {
Ok(Descriptor::P)
}
} else {
Ok(Descriptor::Unknown)
}
}
}
impl CipNode {
pub(crate) const EXPANDED: i32 = 0x1;
pub(crate) const RING_DUPLICATE: i32 = 0x2;
pub(crate) const BOND_DUPLICATE: i32 = 0x4;
pub(crate) const DUPLICATE: i32 = Self::RING_DUPLICATE | Self::BOND_DUPLICATE;
pub(crate) const IMPL_HYDROGEN: i32 = 0x8;
pub(crate) const DUPLICATE_OR_H: i32 =
Self::RING_DUPLICATE | Self::BOND_DUPLICATE | Self::IMPL_HYDROGEN;
pub(crate) const NO_ATOM_INDEX: usize = usize::MAX;
pub(crate) fn new(
digraph: usize,
visit: Vec<i8>,
atom_idx: Option<usize>,
frac: RationalI32,
dist: i32,
flags: i32,
mol: &CipMol<'_>,
) -> Result<Self, CipLabelerError> {
let mut flags = flags;
let (edge_capacity, atomic_mass) = if flags & Self::DUPLICATE != 0 {
(4, 0.0)
} else {
let atomic_number = atom_idx
.map(|idx| mol.atom(idx).map(Atom::atomic_number))
.transpose()?
.unwrap_or(1);
let isotope = atom_idx
.map(|idx| mol.atom(idx).map(Atom::isotope))
.transpose()?
.flatten();
(0, rdkit_atomic_mass(atomic_number, isotope))
};
if visit.is_empty() || flags & Self::DUPLICATE != 0 {
flags |= Self::EXPANDED;
}
Ok(Self {
digraph,
atom_idx,
distance: dist,
atomic_num_fraction: frac,
atomic_mass,
aux: Descriptor::None,
flags,
edges: Vec::with_capacity(edge_capacity),
visit,
})
}
fn new_terminal_child(
&self,
idx: Option<usize>,
atom_idx: Option<usize>,
flags: i32,
mol: &mut CipMol<'_>,
) -> Result<Self, CipLabelerError> {
let new_dist = if flags & Self::DUPLICATE != 0 {
let idx = idx.ok_or(CipLabelerError::UnsupportedDependency {
dependency: "Node::newTerminalChild duplicate without atom index",
})?;
i32::from(self.visit[idx])
} else {
self.distance + 1
};
let new_visit = Vec::new();
if flags & Self::BOND_DUPLICATE != 0 {
let current_atom = self
.atom_idx
.ok_or(CipLabelerError::UnsupportedDependency {
dependency: "Node::newTerminalChild bond duplicate from null atom",
})?;
let frac = mol.get_fractional_atomic_num(current_atom)?;
if frac.denominator > 1 {
return Self::new(
self.digraph,
new_visit,
atom_idx,
frac,
new_dist,
flags,
mol,
);
}
}
let atomic_num = atom_idx
.map(|idx| mol.atom(idx).map(Atom::atomic_number))
.transpose()?
.unwrap_or(1);
Self::new(
self.digraph,
new_visit,
atom_idx,
RationalI32::new(i32::from(atomic_num), 1),
new_dist,
flags,
mol,
)
}
pub(crate) fn get_digraph(&self) -> usize {
self.digraph
}
pub(crate) fn atom_idx(&self) -> Option<usize> {
self.atom_idx
}
pub(crate) fn get_atom_idx(&self) -> usize {
if self.is_set(Self::IMPL_HYDROGEN) {
Self::NO_ATOM_INDEX
} else {
self.atom_idx
.expect("RDKit Node::getAtomIdx requires non-null atom")
}
}
pub(crate) fn get_distance(&self) -> i32 {
self.distance
}
pub(crate) fn get_atomic_num_fraction(&self) -> RationalI32 {
self.atomic_num_fraction
}
pub(crate) fn get_atomic_num(&self, mol: &CipMol<'_>) -> Result<u8, CipLabelerError> {
self.atom_idx
.map(|idx| mol.atom(idx).map(Atom::atomic_number))
.transpose()
.map(|atomic_number| atomic_number.unwrap_or(1))
}
pub(crate) fn get_mass_num(&self, mol: &CipMol<'_>) -> Result<u16, CipLabelerError> {
if self.atom_idx.is_none() || self.is_duplicate() {
return Ok(0);
}
Ok(mol
.atom(self.atom_idx.expect("checked"))?
.isotope()
.unwrap_or(0))
}
pub(crate) fn get_atomic_mass(&self) -> f64 {
self.atomic_mass
}
pub(crate) fn get_aux(&self) -> Descriptor {
self.aux
}
pub(crate) fn is_set(&self, mask: i32) -> bool {
mask & self.flags != 0
}
pub(crate) fn is_duplicate(&self) -> bool {
self.flags & Self::DUPLICATE != 0
}
pub(crate) fn is_duplicate_or_h(&self) -> bool {
self.flags & Self::DUPLICATE_OR_H != 0
}
pub(crate) fn is_terminal(&self) -> bool {
self.visit.is_empty() || (self.is_expanded() && self.edges.len() == 1)
}
pub(crate) fn is_expanded(&self) -> bool {
self.flags & Self::EXPANDED != 0
}
pub(crate) fn is_visited(&self, idx: usize) -> bool {
self.visit[idx] != 0
}
pub(crate) fn new_child(
&self,
idx: usize,
atom_idx: Option<usize>,
mol: &CipMol<'_>,
) -> Result<Self, CipLabelerError> {
let mut new_visit = self.visit.clone();
new_visit[idx] = (self.distance + 1) as i8;
let atomic_num = atom_idx
.map(|idx| mol.atom(idx).map(Atom::atomic_number))
.transpose()?
.unwrap_or(1);
Self::new(
self.digraph,
new_visit,
atom_idx,
RationalI32::new(i32::from(atomic_num), 1),
self.distance + 1,
0,
mol,
)
}
pub(crate) fn new_bond_duplicate_child(
&self,
idx: usize,
atom_idx: Option<usize>,
mol: &mut CipMol<'_>,
) -> Result<Self, CipLabelerError> {
self.new_terminal_child(Some(idx), atom_idx, Self::BOND_DUPLICATE, mol)
}
pub(crate) fn new_ring_duplicate_child(
&self,
idx: usize,
atom_idx: Option<usize>,
mol: &mut CipMol<'_>,
) -> Result<Self, CipLabelerError> {
self.new_terminal_child(Some(idx), atom_idx, Self::RING_DUPLICATE, mol)
}
pub(crate) fn new_implicit_hydrogen_child(
&self,
mol: &mut CipMol<'_>,
) -> Result<Self, CipLabelerError> {
self.new_terminal_child(None, None, Self::IMPL_HYDROGEN, mol)
}
pub(crate) fn add(&mut self, edge: CipEdgeId) {
self.edges.push(edge);
}
pub(crate) fn set_aux(&mut self, desc: Descriptor) {
self.aux = desc;
}
pub(crate) fn get_edges(&self) -> Result<&[CipEdgeId], CipLabelerError> {
if !self.is_expanded() {
return Err(CipLabelerError::UnsupportedDependency {
dependency: "use CipDigraph::node_edges for Node::getEdges lazy expansion",
});
}
Ok(&self.edges)
}
pub(crate) fn get_edges_for_atom(
&self,
end_atom_idx: Option<usize>,
nodes: &[CipNode],
edges: &[CipEdge],
) -> Result<Vec<CipEdgeId>, CipLabelerError> {
let mut result = Vec::new();
for edge_id in self.get_edges()? {
let edge = &edges[edge_id.index()];
if nodes[edge.end.index()].is_duplicate() {
continue;
}
if end_atom_idx == nodes[edge.beg.index()].atom_idx
|| end_atom_idx == nodes[edge.end.index()].atom_idx
{
result.push(*edge_id);
}
}
Ok(result)
}
pub(crate) fn get_non_terminal_out_edges(
&self,
self_id: CipNodeId,
nodes: &[CipNode],
edges: &[CipEdge],
) -> Result<Vec<CipEdgeId>, CipLabelerError> {
let mut result = Vec::new();
for edge_id in self.get_edges()? {
let edge = &edges[edge_id.index()];
if edge.beg == self_id && !nodes[edge.end.index()].is_terminal() {
result.push(*edge_id);
}
}
Ok(result)
}
}
impl<'a> CipDigraph<'a> {
const MAX_NODE_COUNT: usize = 100_000;
const MAX_NODE_DIST: i32 = 0;
pub(crate) fn new(
molecule: &'a Molecule,
atom_idx: usize,
atropisomer_mode: bool,
) -> Result<Self, CipLabelerError> {
let mol = CipMol::new(molecule);
mol.atom(atom_idx)?;
let mut digraph = Self {
mol,
origin: CipNodeId::new(0),
root: CipNodeId::new(0),
rule6_ref: None,
atropisomer_mode,
nodes: Vec::new(),
edges: Vec::new(),
};
let mut visit = vec![0_i8; digraph.mol.get_num_atoms()];
visit[atom_idx] = 1;
let atomic_num = digraph.mol.atom(atom_idx)?.atomic_number();
let root = digraph.add_node(
visit,
Some(atom_idx),
RationalI32::new(i32::from(atomic_num), 1),
1,
0,
)?;
digraph.root = root;
digraph.origin = root;
Ok(digraph)
}
pub(crate) fn mol(&self) -> &CipMol<'a> {
&self.mol
}
pub(crate) fn get_original_root(&self) -> CipNodeId {
self.origin
}
pub(crate) fn get_current_root(&self) -> CipNodeId {
self.root
}
pub(crate) fn get_num_nodes(&self) -> usize {
self.nodes.len()
}
pub(crate) fn node(&self, node: CipNodeId) -> &CipNode {
&self.nodes[node.index()]
}
pub(crate) fn edge(&self, edge: CipEdgeId) -> &CipEdge {
&self.edges[edge.index()]
}
fn add_node(
&mut self,
visit: Vec<i8>,
atom_idx: Option<usize>,
frac: RationalI32,
dist: i32,
flags: i32,
) -> Result<CipNodeId, CipLabelerError> {
let node = CipNode::new(0, visit, atom_idx, frac, dist, flags, &self.mol)?;
let id = CipNodeId::new(self.nodes.len());
self.nodes.push(node);
Ok(id)
}
fn add_edge(&mut self, beg: CipNodeId, bond_idx: Option<usize>, end: CipNodeId) {
let edge_id = CipEdgeId::new(self.edges.len());
self.edges.push(CipEdge::new(beg, end, bond_idx));
self.nodes[beg.index()].add(edge_id);
self.nodes[end.index()].add(edge_id);
}
pub(crate) fn node_edges(
&mut self,
node: CipNodeId,
) -> Result<Vec<CipEdgeId>, CipLabelerError> {
if !self.nodes[node.index()].is_expanded() {
self.nodes[node.index()].flags |= CipNode::EXPANDED;
self.expand(node)?;
}
Ok(self.nodes[node.index()].edges.clone())
}
pub(crate) fn node_edges_for_atom(
&mut self,
node: CipNodeId,
end_atom_idx: Option<usize>,
) -> Result<Vec<CipEdgeId>, CipLabelerError> {
let edge_ids = self.node_edges(node)?;
let mut result = Vec::new();
for edge_id in edge_ids {
let edge = &self.edges[edge_id.index()];
if self.nodes[edge.get_end().index()].is_duplicate() {
continue;
}
if end_atom_idx == self.nodes[edge.get_beg().index()].atom_idx
|| end_atom_idx == self.nodes[edge.get_end().index()].atom_idx
{
result.push(edge_id);
}
}
Ok(result)
}
pub(crate) fn non_terminal_out_edges(
&mut self,
node: CipNodeId,
) -> Result<Vec<CipEdgeId>, CipLabelerError> {
let edge_ids = self.node_edges(node)?;
let mut result = Vec::new();
for edge_id in edge_ids {
let edge = &self.edges[edge_id.index()];
if edge.is_beg(node) && !self.nodes[edge.get_end().index()].is_terminal() {
result.push(edge_id);
}
}
Ok(result)
}
pub(crate) fn get_nodes(&mut self, atom_idx: usize) -> Result<Vec<CipNodeId>, CipLabelerError> {
self.mol.atom(atom_idx)?;
let mut result = Vec::new();
let mut queue = vec![self.get_current_root()];
let mut i = 0_usize;
while i < queue.len() {
let node = queue[i];
if self.nodes[node.index()].atom_idx == Some(atom_idx) {
result.push(node);
}
for edge_id in self.node_edges(node)? {
let edge = &self.edges[edge_id.index()];
if !edge.is_beg(node) {
continue;
}
queue.push(edge.get_end());
}
i += 1;
}
Ok(result)
}
pub(crate) fn get_rule6_ref(&self) -> Option<usize> {
self.rule6_ref
}
pub(crate) fn set_rule6_ref(&mut self, atom_idx: Option<usize>) -> Result<(), CipLabelerError> {
if let Some(atom_idx) = atom_idx {
self.mol.atom(atom_idx)?;
}
self.rule6_ref = atom_idx;
Ok(())
}
pub(crate) fn change_root(&mut self, new_root: CipNodeId) -> Result<(), CipLabelerError> {
let mut to_flip = Vec::new();
let mut queue = vec![new_root];
let mut i = 0_usize;
while i < queue.len() {
let node = queue[i];
for edge_id in self.node_edges(node)? {
let edge = &self.edges[edge_id.index()];
if edge.is_end(node) {
to_flip.push(edge_id);
queue.push(edge.get_beg());
}
}
i += 1;
}
for edge_id in to_flip {
self.edges[edge_id.index()].flip();
}
self.root = new_root;
Ok(())
}
fn expand(&mut self, beg: CipNodeId) -> Result<(), CipLabelerError> {
let atom_idx =
self.nodes[beg.index()]
.atom_idx
.ok_or(CipLabelerError::UnsupportedDependency {
dependency: "Digraph::expand on null atom node",
})?;
let prev = self.nodes[beg.index()].edges.first().and_then(|edge_id| {
let edge = &self.edges[edge_id.index()];
(!edge.is_beg(beg)).then_some(edge.get_bond_idx()).flatten()
});
if Self::MAX_NODE_DIST > 0 && self.nodes[beg.index()].get_distance() > Self::MAX_NODE_DIST {
return Ok(());
}
if Self::MAX_NODE_COUNT > 0 && self.nodes.len() >= Self::MAX_NODE_COUNT {
return Err(CipLabelerError::TooManyNodes {
limit: Self::MAX_NODE_COUNT,
});
}
for bond_idx in self.mol.bond_indices_for_atom(atom_idx)? {
let nbr_idx = self.mol.other_atom_idx(bond_idx, atom_idx)?;
let bord = self.mol.get_bond_order(bond_idx)?;
let virtual_nodes = bord - 1;
if !self.nodes[beg.index()].is_visited(nbr_idx) {
let end_node =
self.nodes[beg.index()].new_child(nbr_idx, Some(nbr_idx), &self.mol)?;
let mut end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
if self.origin != beg || self.atropisomer_mode {
let atom_formal_charge = self.mol.atom(atom_idx)?.formal_charge();
if atom_formal_charge < 0
&& self.mol.get_fractional_atomic_num(atom_idx)?.denominator > 1
{
let end_node = self.nodes[beg.index()].new_bond_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
} else {
for _ in 0..virtual_nodes {
let end_node = self.nodes[beg.index()].new_bond_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
}
}
}
} else if Some(bond_idx) == prev {
if self.nodes[self.origin.index()].atom_idx != Some(nbr_idx)
|| self.atropisomer_mode
{
for _ in 0..virtual_nodes {
let end_node = self.nodes[beg.index()].new_bond_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
let end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
}
}
} else {
let end_node = self.nodes[beg.index()].new_ring_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
let mut end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
let atom_formal_charge = self.mol.atom(atom_idx)?.formal_charge();
if atom_formal_charge < 0
&& self.mol.get_fractional_atomic_num(atom_idx)?.denominator > 1
{
let end_node = self.nodes[beg.index()].new_bond_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
} else {
for _ in 0..virtual_nodes {
let end_node = self.nodes[beg.index()].new_bond_duplicate_child(
nbr_idx,
Some(nbr_idx),
&mut self.mol,
)?;
end = self.add_existing_node(end_node);
self.add_edge(beg, Some(bond_idx), end);
}
}
}
}
let hcnt = self.mol.total_num_hs(atom_idx)?;
for _ in 0..hcnt {
let end_node = self.nodes[beg.index()].new_implicit_hydrogen_child(&mut self.mol)?;
let end = self.add_existing_node(end_node);
self.add_edge(beg, None, end);
}
Ok(())
}
fn add_existing_node(&mut self, node: CipNode) -> CipNodeId {
let id = CipNodeId::new(self.nodes.len());
self.nodes.push(node);
id
}
}
impl<'a> CipMol<'a> {
pub(crate) fn new(molecule: &'a Molecule) -> Self {
Self {
molecule,
rings: None,
kekulized_bond_orders: None,
fractional_atomic_numbers: None,
valence: None,
}
}
pub(crate) fn get_fractional_atomic_num(
&mut self,
atom_idx: usize,
) -> Result<RationalI32, CipLabelerError> {
self.atom(atom_idx)?;
if self.fractional_atomic_numbers.is_none() {
self.fractional_atomic_numbers = Some(calc_frac_atom_nums(self)?);
}
Ok(self
.fractional_atomic_numbers
.as_ref()
.expect("initialized")[atom_idx])
}
pub(crate) fn get_num_atoms(&self) -> usize {
self.molecule.num_atoms()
}
pub(crate) fn get_num_bonds(&self) -> usize {
self.molecule.num_bonds()
}
pub(crate) fn atom(&self, atom_idx: usize) -> Result<&'a Atom, CipLabelerError> {
self.molecule
.atoms()
.get(atom_idx)
.ok_or(CipLabelerError::AtomIndexOutOfRange {
index: atom_idx,
atom_count: self.molecule.num_atoms(),
})
}
pub(crate) fn atoms(&self) -> &'a [Atom] {
self.molecule.atoms()
}
pub(crate) fn bond(&self, bond_idx: usize) -> Result<&'a Bond, CipLabelerError> {
self.molecule
.bonds()
.get(bond_idx)
.ok_or(CipLabelerError::BondIndexOutOfRange {
index: bond_idx,
bond_count: self.molecule.num_bonds(),
})
}
pub(crate) fn bond_indices_for_atom(
&self,
atom_idx: usize,
) -> Result<Vec<usize>, CipLabelerError> {
self.atom(atom_idx)?;
Ok(self
.molecule
.topology_block()
.adjacency
.neighbors_of(atom_idx)
.iter()
.map(|neighbor| neighbor.bond.index())
.collect())
}
pub(crate) fn neighbor_indices(&self, atom_idx: usize) -> Result<Vec<usize>, CipLabelerError> {
self.atom(atom_idx)?;
Ok(self
.molecule
.topology_block()
.adjacency
.neighbors_of(atom_idx)
.iter()
.map(|neighbor| neighbor.atom_index)
.collect())
}
pub(crate) fn is_in_ring(&mut self, bond_idx: usize) -> Result<bool, CipLabelerError> {
self.bond(bond_idx)?;
if self.rings.is_none() {
self.rings = Some(crate::rings::fast_find_rings(self.molecule)?);
}
Ok(self
.rings
.as_ref()
.expect("initialized")
.num_bond_rings(BondId::new(bond_idx))
!= 0)
}
pub(crate) fn get_bond_order(&mut self, bond_idx: usize) -> Result<i32, CipLabelerError> {
self.bond(bond_idx)?;
if self.kekulized_bond_orders.is_none() {
let mut orders = self
.molecule
.bonds()
.iter()
.map(Bond::order)
.collect::<Vec<_>>();
if let Ok(assignment) = crate::kekulize::kekulize_assignment(
self.molecule,
self.rings.as_ref(),
true,
true,
100,
) {
for (idx, order) in orders.iter_mut().enumerate() {
if let Some(kekulized) = assignment.bond_order(BondId::new(idx)) {
*order = kekulized;
}
}
}
self.kekulized_bond_orders = Some(orders);
}
match self.kekulized_bond_orders.as_ref().expect("initialized")[bond_idx] {
BondOrder::Zero
| BondOrder::Hydrogen
| BondOrder::Dative
| BondOrder::DativeLeft
| BondOrder::DativeRight => Ok(0),
BondOrder::Single => Ok(1),
BondOrder::Aromatic => Ok(1),
BondOrder::Double => Ok(2),
BondOrder::Triple => Ok(3),
BondOrder::Quadruple => Ok(4),
BondOrder::Quintuple => Ok(5),
BondOrder::Hextuple => Ok(6),
order => Err(CipLabelerError::NonIntegerBondOrder { order }),
}
}
fn total_num_hs(&mut self, atom_idx: usize) -> Result<i32, CipLabelerError> {
let explicit = i32::from(self.atom(atom_idx)?.explicit_hydrogens());
if self.valence.is_none() {
self.valence = Some(crate::valence::assign_valence(
self.molecule,
ValenceModel::RdkitLike,
)?);
}
let implicit = self
.valence
.as_ref()
.and_then(|valence| valence.implicit_hydrogens.get(atom_idx))
.copied()
.unwrap_or(0)
.max(0);
Ok(explicit + implicit)
}
fn other_atom_idx(&self, bond_idx: usize, atom_idx: usize) -> Result<usize, CipLabelerError> {
let bond = self.bond(bond_idx)?;
if bond.begin().index() == atom_idx {
Ok(bond.end().index())
} else if bond.end().index() == atom_idx {
Ok(bond.begin().index())
} else {
Err(CipLabelerError::BondNotIncident {
bond: bond_idx,
atom: atom_idx,
})
}
}
}
fn seed_types(types: &mut [MancudeType], mol: &mut CipMol<'_>) -> Result<bool, CipLabelerError> {
let mut result = false;
for atom_idx in 0..mol.get_num_atoms() {
let mut btypes = mol.total_num_hs(atom_idx)?;
let mut ring = false;
for bond_idx in mol.bond_indices_for_atom(atom_idx)? {
match mol.get_bond_order(bond_idx)? {
1 => btypes += 0x00000001,
2 => btypes += 0x00000100,
_ => btypes += 0x01000000,
}
if mol.is_in_ring(bond_idx)? {
ring = true;
}
}
if !ring {
continue;
}
let atom = mol.atom(atom_idx)?;
let q = atom.formal_charge();
match atom.atomic_number() {
6 | 14 | 32 => {
if q == 0 && btypes == 0x0102 {
types[atom_idx] = MancudeType::Cv4D3;
} else if q == -1 && btypes == 0x0003 {
types[atom_idx] = MancudeType::Cv3D3Minus;
result = true;
}
}
7 | 15 | 33 => {
if q == 0 && btypes == 0x0101 {
types[atom_idx] = MancudeType::Nv3D2;
result = true;
} else if q == -1 && btypes == 0x0002 {
types[atom_idx] = MancudeType::Nv2D2Minus;
result = true;
} else if q == 1 && btypes == 0x0102 {
types[atom_idx] = MancudeType::Nv4D3Plus;
result = true;
}
}
8 => {
if q == 1 && btypes == 0x0101 {
types[atom_idx] = MancudeType::Ov3D2Plus;
result = true;
}
}
_ => {}
}
}
Ok(result)
}
fn relax_types(types: &mut [MancudeType], mol: &CipMol<'_>) -> Result<(), CipLabelerError> {
let mut queue = VecDeque::new();
let mut counts = vec![0_i32; mol.get_num_atoms()];
for atom_idx in 0..mol.get_num_atoms() {
for nbr_idx in mol.neighbor_indices(atom_idx)? {
if types[nbr_idx] != MancudeType::Other {
counts[atom_idx] += 1;
}
}
if counts[atom_idx] == 1 {
queue.push_back(atom_idx);
}
}
while let Some(atom_idx) = queue.pop_front() {
if types[atom_idx] == MancudeType::Other {
continue;
}
types[atom_idx] = MancudeType::Other;
for nbr_idx in mol.neighbor_indices(atom_idx)? {
counts[nbr_idx] -= 1;
if counts[nbr_idx] == 1 {
queue.push_back(nbr_idx);
}
}
}
Ok(())
}
fn visit_part(
parts: &mut [i32],
types: &[MancudeType],
part: i32,
mut atom_idx: usize,
mol: &mut CipMol<'_>,
) -> Result<(), CipLabelerError> {
loop {
let mut next = None;
for bond_idx in mol.bond_indices_for_atom(atom_idx)? {
if !mol.is_in_ring(bond_idx)? {
continue;
}
let nbr_idx = mol.other_atom_idx(bond_idx, atom_idx)?;
if parts[nbr_idx] == 0 && types[nbr_idx] != MancudeType::Other {
parts[nbr_idx] = part;
if next.is_some() {
visit_part(parts, types, part, nbr_idx, mol)?;
} else {
next = Some(nbr_idx);
}
}
}
if let Some(next_idx) = next {
atom_idx = next_idx;
} else {
break;
}
}
Ok(())
}
fn visit_parts(
parts: &mut [i32],
types: &[MancudeType],
mol: &mut CipMol<'_>,
) -> Result<i32, CipLabelerError> {
let mut numparts = 0_i32;
for atom_idx in 0..mol.get_num_atoms() {
if parts[atom_idx] == 0 && types[atom_idx] != MancudeType::Other {
numparts += 1;
parts[atom_idx] = numparts;
visit_part(parts, types, numparts, atom_idx, mol)?;
}
}
Ok(numparts)
}
fn calc_frac_atom_nums(mol: &mut CipMol<'_>) -> Result<Vec<RationalI32>, CipLabelerError> {
let num_atoms = mol.get_num_atoms();
let mut fractions = Vec::with_capacity(num_atoms);
for atom_idx in 0..num_atoms {
fractions.push(RationalI32::new(
i32::from(mol.atom(atom_idx)?.atomic_number()),
1,
));
}
let mut types = vec![MancudeType::Other; num_atoms];
if seed_types(&mut types, mol)? {
relax_types(&mut types, mol)?;
let mut parts = vec![0_i32; num_atoms];
let numparts = visit_parts(&mut parts, &types, mol)?;
let mut resparts = vec![0_i32; usize::try_from(numparts).unwrap_or(0)];
let mut numres = 0_usize;
if numparts > 0 {
for i in 0..num_atoms {
if parts[i] == 0 {
continue;
}
if matches!(types[i], MancudeType::Cv3D3Minus | MancudeType::Nv2D2Minus) {
let mut j = 0_usize;
while j < numres {
if resparts[j] == parts[i] {
break;
}
j += 1;
}
if j >= numres {
resparts[numres] = parts[i];
numres += 1;
}
}
let mut numerator = 0_i32;
let mut denominator = 0_i32;
for nbr_idx in mol.neighbor_indices(i)? {
if parts[nbr_idx] == parts[i] {
numerator += i32::from(mol.atom(nbr_idx)?.atomic_number());
denominator += 1;
}
}
if denominator == 0 {
fractions[i].assign(0, 1);
} else {
fractions[i].assign(numerator, denominator);
}
}
}
if numres > 0 {
for &part in resparts.iter().take(numres) {
let mut numerator = 0_i32;
let mut denominator = 0_i32;
for i in 0..num_atoms {
if parts[i] != part {
continue;
}
if denominator == 0 {
fractions[i].assign(0, 1);
} else {
fractions[i].assign(numerator, denominator);
}
denominator += 1;
for bond_idx in mol.bond_indices_for_atom(i)? {
let nbr_idx = mol.other_atom_idx(bond_idx, i)?;
let bord = mol.get_bond_order(bond_idx)?;
if bord > 1 && parts[nbr_idx] == part {
numerator += (bord - 1) * i32::from(mol.atom(nbr_idx)?.atomic_number());
}
}
}
}
}
}
Ok(fractions)
}
#[cfg(test)]
mod tests {
use super::{
CipAtropisomerBond, CipConfiguration, CipDigraph, CipEdge, CipEdgeId, CipLabelerContext,
CipLabelerError, CipMol, CipNode, CipNodeId, CipPairList, CipPriority, CipRule1a,
CipRule1b, CipRule2, CipRule3, CipRule4a, CipRule4b, CipRule4c, CipRule5New, CipRule6,
CipRules, CipSequenceRule, CipSort, CipSp2Bond, CipTetrahedral, Descriptor, RationalI32,
assign_cip_labels, assign_cip_labels_for_indices, cip_all_rules, descriptor_to_string,
};
use crate::{
AtomSpec, BondOrder, BondSpec, BondStereo, ChiralTag, Element, Molecule, MoleculeBuilder,
};
#[test]
fn ciplabeler_descriptor_to_string_matches_rdkit() {
let expected = [
(Descriptor::None, "NONE"),
(Descriptor::Unknown, "UNKNOWN"),
(Descriptor::ns, "ns"),
(Descriptor::R, "R"),
(Descriptor::S, "S"),
(Descriptor::r, "r"),
(Descriptor::s, "s"),
(Descriptor::seqTrans, "e"),
(Descriptor::seqCis, "z"),
(Descriptor::E, "E"),
(Descriptor::Z, "Z"),
(Descriptor::M, "M"),
(Descriptor::P, "P"),
(Descriptor::m, "m"),
(Descriptor::p, "p"),
(Descriptor::SP_4, "SP_4"),
(Descriptor::TBPY_5, "TBPY_5"),
(Descriptor::OC_6, "OC_6"),
];
for (desc, label) in expected {
assert_eq!(descriptor_to_string(desc), label);
}
}
#[test]
fn ciplabeler_descriptor_order_matches_rdkit_enum_order() {
assert_eq!(
Descriptor::ALL_IN_RDKIT_ORDER,
[
Descriptor::None,
Descriptor::Unknown,
Descriptor::ns,
Descriptor::R,
Descriptor::S,
Descriptor::r,
Descriptor::s,
Descriptor::seqTrans,
Descriptor::seqCis,
Descriptor::E,
Descriptor::Z,
Descriptor::M,
Descriptor::P,
Descriptor::m,
Descriptor::p,
Descriptor::SP_4,
Descriptor::TBPY_5,
Descriptor::OC_6,
]
);
for window in Descriptor::ALL_IN_RDKIT_ORDER.windows(2) {
assert!(window[0] < window[1]);
}
}
#[test]
fn ciplabeler_cipmol_accessors_match_rdkit_molecule_view() {
let molecule = Molecule::from_smiles("CCO").unwrap();
let cipmol = CipMol::new(&molecule);
assert_eq!(cipmol.get_num_atoms(), 3);
assert_eq!(cipmol.get_num_bonds(), 2);
assert_eq!(cipmol.atom(0).unwrap().atomic_number(), 6);
assert_eq!(cipmol.atom(2).unwrap().atomic_number(), 8);
assert_eq!(cipmol.bond(0).unwrap().order(), BondOrder::Single);
assert_eq!(cipmol.neighbor_indices(1).unwrap(), vec![0, 2]);
assert_eq!(cipmol.bond_indices_for_atom(1).unwrap(), vec![0, 1]);
}
#[test]
fn ciplabeler_cipmol_ring_membership_uses_fast_ring_info() {
let molecule = Molecule::from_smiles("C1CCCCC1").unwrap();
let mut cipmol = CipMol::new(&molecule);
assert!(cipmol.is_in_ring(0).unwrap());
assert!(cipmol.is_in_ring(5).unwrap());
let chain = Molecule::from_smiles("CCCC").unwrap();
let mut chain_cipmol = CipMol::new(&chain);
assert!(!chain_cipmol.is_in_ring(0).unwrap());
}
#[test]
fn ciplabeler_cipmol_bond_order_uses_kekulized_aromatic_copy() {
let molecule = Molecule::from_smiles("c1ccccc1").unwrap();
let mut cipmol = CipMol::new(&molecule);
let mut orders = (0..molecule.num_bonds())
.map(|idx| cipmol.get_bond_order(idx).unwrap())
.collect::<Vec<_>>();
orders.sort_unstable();
assert_eq!(orders, vec![1, 1, 1, 2, 2, 2]);
}
#[test]
fn ciplabeler_cipmol_bond_order_maps_zero_hydrogen_and_dative_to_zero() {
let mut builder = MoleculeBuilder::new();
let a = builder.add_atom(AtomSpec::new(Element::N));
let b = builder.add_atom(AtomSpec::new(Element::B));
let c = builder.add_atom(AtomSpec::new(Element::H));
let d = builder.add_atom(AtomSpec::new(Element::H));
builder
.add_bond(BondSpec::new(a, b, BondOrder::Dative))
.unwrap();
builder
.add_bond(BondSpec::new(b, c, BondOrder::Zero))
.unwrap();
builder
.add_bond(BondSpec::new(c, d, BondOrder::Hydrogen))
.unwrap();
let molecule = builder.build().unwrap();
let mut cipmol = CipMol::new(&molecule);
assert_eq!(cipmol.get_bond_order(0).unwrap(), 0);
assert_eq!(cipmol.get_bond_order(1).unwrap(), 0);
assert_eq!(cipmol.get_bond_order(2).unwrap(), 0);
}
#[test]
fn ciplabeler_cipmol_fractional_atomic_numbers_handle_mancude_nitrogen() {
let molecule = Molecule::from_smiles("c1ccncc1").unwrap();
let mut cipmol = CipMol::new(&molecule);
let fractions = (0..molecule.num_atoms())
.map(|idx| cipmol.get_fractional_atomic_num(idx).unwrap().tuple())
.collect::<Vec<_>>();
assert_eq!(
fractions,
vec![(6, 1), (6, 1), (13, 2), (6, 1), (13, 2), (6, 1)]
);
}
#[test]
fn ciplabeler_cipmol_fractional_atomic_numbers_recompute_negative_charge_resonance() {
let molecule = Molecule::from_smiles("C1=C[CH-]C=C1").unwrap();
let mut cipmol = CipMol::new(&molecule);
let fractions = (0..molecule.num_atoms())
.map(|idx| cipmol.get_fractional_atomic_num(idx).unwrap().tuple())
.collect::<Vec<_>>();
assert_eq!(fractions, vec![(0, 1), (6, 1), (6, 1), (4, 1), (9, 2)]);
}
#[test]
fn ciplabeler_node_constructor_flags_mass_and_visits_match_rdkit() {
let molecule = Molecule::from_smiles("[13CH4]").unwrap();
let cipmol = CipMol::new(&molecule);
let node =
CipNode::new(11, vec![1], Some(0), RationalI32::new(6, 1), 1, 0, &cipmol).unwrap();
assert_eq!(node.get_digraph(), 11);
assert_eq!(node.atom_idx(), Some(0));
assert_eq!(node.get_atom_idx(), 0);
assert_eq!(node.get_distance(), 1);
assert_eq!(node.get_atomic_num_fraction().tuple(), (6, 1));
assert_eq!(node.get_atomic_num(&cipmol).unwrap(), 6);
assert_eq!(node.get_mass_num(&cipmol).unwrap(), 13);
assert!(node.get_atomic_mass() > 13.0);
assert_eq!(node.get_aux(), Descriptor::None);
assert!(!node.is_duplicate());
assert!(!node.is_duplicate_or_h());
assert!(!node.is_expanded());
assert!(node.is_visited(0));
assert!(!node.is_terminal());
let duplicate = CipNode::new(
11,
vec![1],
Some(0),
RationalI32::new(6, 1),
1,
CipNode::BOND_DUPLICATE,
&cipmol,
)
.unwrap();
assert!(duplicate.is_duplicate());
assert!(duplicate.is_duplicate_or_h());
assert!(duplicate.is_expanded());
assert_eq!(duplicate.get_mass_num(&cipmol).unwrap(), 0);
assert_eq!(duplicate.get_atomic_mass(), 0.0);
}
#[test]
fn ciplabeler_node_child_creation_and_aux_storage_match_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut cipmol = CipMol::new(&molecule);
let root = CipNode::new(
3,
vec![1, 0],
Some(0),
RationalI32::new(6, 1),
1,
0,
&cipmol,
)
.unwrap();
let child = root.new_child(1, Some(1), &cipmol).unwrap();
assert_eq!(child.atom_idx(), Some(1));
assert_eq!(child.get_distance(), 2);
assert!(child.is_visited(0));
assert!(child.is_visited(1));
assert!(!child.is_expanded());
let ring_duplicate = root
.new_ring_duplicate_child(0, Some(0), &mut cipmol)
.unwrap();
assert!(ring_duplicate.is_duplicate());
assert!(ring_duplicate.is_expanded());
assert_eq!(ring_duplicate.get_distance(), 1);
assert_eq!(ring_duplicate.get_atomic_num_fraction().tuple(), (6, 1));
let implicit_h = root.new_implicit_hydrogen_child(&mut cipmol).unwrap();
assert_eq!(implicit_h.atom_idx(), None);
assert_eq!(implicit_h.get_atom_idx(), CipNode::NO_ATOM_INDEX);
assert_eq!(implicit_h.get_atomic_num(&cipmol).unwrap(), 1);
assert_eq!(implicit_h.get_mass_num(&cipmol).unwrap(), 0);
assert!(implicit_h.is_duplicate_or_h());
assert!(implicit_h.is_expanded());
let mut aux_node = root.clone();
aux_node.set_aux(Descriptor::R);
assert_eq!(aux_node.get_aux(), Descriptor::R);
}
#[test]
fn ciplabeler_edge_endpoints_aux_and_flip_match_rdkit() {
let mut edge = CipEdge::new(CipNodeId::new(0), CipNodeId::new(1), Some(7));
assert_eq!(edge.get_beg(), CipNodeId::new(0));
assert_eq!(edge.get_end(), CipNodeId::new(1));
assert_eq!(edge.get_bond_idx(), Some(7));
assert!(edge.is_beg(CipNodeId::new(0)));
assert!(edge.is_end(CipNodeId::new(1)));
assert_eq!(
edge.get_other(CipEdgeId::new(4), CipNodeId::new(0))
.unwrap(),
CipNodeId::new(1)
);
assert_eq!(
edge.get_other(CipEdgeId::new(4), CipNodeId::new(1))
.unwrap(),
CipNodeId::new(0)
);
assert!(matches!(
edge.get_other(CipEdgeId::new(4), CipNodeId::new(2)),
Err(CipLabelerError::EdgeEndpointMismatch { edge: 4, node: 2 })
));
edge.set_aux(Descriptor::seqTrans);
assert_eq!(edge.get_aux(), Descriptor::seqTrans);
edge.flip();
assert_eq!(edge.get_beg(), CipNodeId::new(1));
assert_eq!(edge.get_end(), CipNodeId::new(0));
}
#[test]
fn ciplabeler_node_raw_edges_fail_closed_until_digraph_lazy_expansion() {
let molecule = Molecule::from_smiles("CC").unwrap();
let cipmol = CipMol::new(&molecule);
let node = CipNode::new(
3,
vec![1, 0],
Some(0),
RationalI32::new(6, 1),
1,
0,
&cipmol,
)
.unwrap();
assert!(matches!(
node.get_edges(),
Err(CipLabelerError::UnsupportedDependency {
dependency: "use CipDigraph::node_edges for Node::getEdges lazy expansion"
})
));
}
#[test]
fn ciplabeler_digraph_root_construction_and_rule6_ref_match_rdkit() {
let molecule = Molecule::from_smiles("CCO").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
assert_eq!(digraph.get_original_root(), CipNodeId::new(0));
assert_eq!(digraph.get_current_root(), CipNodeId::new(0));
assert_eq!(digraph.get_num_nodes(), 1);
assert_eq!(digraph.node(digraph.get_current_root()).atom_idx(), Some(1));
assert_eq!(digraph.node(digraph.get_current_root()).get_distance(), 1);
assert_eq!(digraph.mol().get_num_atoms(), 3);
assert_eq!(digraph.get_rule6_ref(), None);
digraph.set_rule6_ref(Some(2)).unwrap();
assert_eq!(digraph.get_rule6_ref(), Some(2));
digraph.set_rule6_ref(None).unwrap();
assert_eq!(digraph.get_rule6_ref(), None);
}
#[test]
fn ciplabeler_digraph_lazy_expansion_creates_explicit_and_implicit_nodes() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
assert_eq!(root_edges.len(), 4);
assert_eq!(digraph.get_num_nodes(), 5);
assert!(digraph.node(root).is_expanded());
let explicit_edges = root_edges
.iter()
.filter(|edge_id| digraph.edge(**edge_id).get_bond_idx().is_some())
.copied()
.collect::<Vec<_>>();
assert_eq!(explicit_edges.len(), 1);
let explicit_end = digraph.edge(explicit_edges[0]).get_end();
assert_eq!(digraph.node(explicit_end).atom_idx(), Some(1));
assert_eq!(digraph.node(explicit_end).get_distance(), 2);
let implicit_count = root_edges
.iter()
.filter(|edge_id| digraph.edge(**edge_id).get_bond_idx().is_none())
.count();
assert_eq!(implicit_count, 3);
}
#[test]
fn ciplabeler_digraph_expansion_creates_bond_and_ring_duplicates() {
let ethene = Molecule::from_smiles("C=C").unwrap();
let mut ethene_graph = CipDigraph::new(ðene, 0, true).unwrap();
let root = ethene_graph.get_current_root();
let root_edges = ethene_graph.node_edges(root).unwrap();
let duplicate_edges = root_edges
.iter()
.filter(|edge_id| {
let end = ethene_graph.edge(**edge_id).get_end();
ethene_graph.node(end).is_duplicate()
})
.count();
assert_eq!(duplicate_edges, 1);
let cyclopropane = Molecule::from_smiles("C1CC1").unwrap();
let mut ring_graph = CipDigraph::new(&cyclopropane, 0, false).unwrap();
let root = ring_graph.get_current_root();
let first_layer = ring_graph.node_edges(root).unwrap();
let non_root_neighbor = first_layer
.iter()
.map(|edge_id| ring_graph.edge(*edge_id).get_end())
.find(|node_id| ring_graph.node(*node_id).atom_idx() == Some(1))
.unwrap();
let second_layer = ring_graph.node_edges(non_root_neighbor).unwrap();
let closing_node = second_layer
.iter()
.map(|edge_id| ring_graph.edge(*edge_id).get_end())
.find(|node_id| ring_graph.node(*node_id).atom_idx() == Some(2))
.unwrap();
let third_layer = ring_graph.node_edges(closing_node).unwrap();
assert!(third_layer.iter().any(|edge_id| {
let end = ring_graph.edge(*edge_id).get_end();
ring_graph.node(end).is_set(CipNode::RING_DUPLICATE)
}));
}
#[test]
fn ciplabeler_digraph_get_nodes_and_change_root_match_rdkit_direction_rules() {
let molecule = Molecule::from_smiles("CCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.map(|edge_id| digraph.edge(*edge_id).get_end())
.find(|node_id| digraph.node(*node_id).atom_idx() == Some(0))
.unwrap();
digraph.node_edges(left).unwrap();
let center_nodes = digraph.get_nodes(1).unwrap();
assert!(center_nodes.iter().any(|node| *node == root));
digraph.change_root(left).unwrap();
assert_eq!(digraph.get_current_root(), left);
let flipped_back_edge = digraph
.node_edges(left)
.unwrap()
.into_iter()
.find(|edge_id| digraph.edge(*edge_id).get_end() == root)
.unwrap();
assert!(digraph.edge(flipped_back_edge).is_beg(left));
}
struct AtomicNumberSortRule;
impl CipSequenceRule for AtomicNumberSortRule {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let a_atomic_num = i32::from(digraph.node(a_end).get_atomic_num(digraph.mol())?);
let b_atomic_num = i32::from(digraph.node(b_end).get_atomic_num(digraph.mol())?);
Ok(match a_atomic_num.cmp(&b_atomic_num) {
std::cmp::Ordering::Less => -1,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => 1,
})
}
}
struct PseudoAsymSortRule;
impl CipSequenceRule for PseudoAsymSortRule {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let a_atomic_num = i32::from(digraph.node(a_end).get_atomic_num(digraph.mol())?);
let b_atomic_num = i32::from(digraph.node(b_end).get_atomic_num(digraph.mol())?);
Ok(match a_atomic_num.cmp(&b_atomic_num) {
std::cmp::Ordering::Less => -2,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => 2,
})
}
}
struct SinglePseudoAsymSortRule;
impl CipSequenceRule for SinglePseudoAsymSortRule {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_end = digraph.edge(a).get_end();
let b_end = digraph.edge(b).get_end();
let a_atomic_num = i32::from(digraph.node(a_end).get_atomic_num(digraph.mol())?);
let b_atomic_num = i32::from(digraph.node(b_end).get_atomic_num(digraph.mol())?);
let magnitude = if (a_atomic_num == 35 && b_atomic_num == 53)
|| (a_atomic_num == 53 && b_atomic_num == 35)
{
2
} else {
1
};
Ok(match a_atomic_num.cmp(&b_atomic_num) {
std::cmp::Ordering::Less => -magnitude,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => magnitude,
})
}
}
struct BondLabelSortRule;
impl CipSequenceRule for BondLabelSortRule {
fn compare(
&self,
digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
a: CipEdgeId,
b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
let a_label = self.get_bond_label(digraph.edge(a));
let b_label = self.get_bond_label(digraph.edge(b));
Ok(match a_label.cmp(&b_label) {
std::cmp::Ordering::Less => -1,
std::cmp::Ordering::Equal => 0,
std::cmp::Ordering::Greater => 1,
})
}
}
struct AlwaysEqualSortRule;
impl CipSequenceRule for AlwaysEqualSortRule {
fn compare(
&self,
_digraph: &mut CipDigraph<'_>,
_context: &mut CipLabelerContext,
_a: CipEdgeId,
_b: CipEdgeId,
) -> Result<i32, CipLabelerError> {
Ok(0)
}
}
fn edge_end_atomic_nums(
digraph: &CipDigraph<'_>,
edges: &[CipEdgeId],
) -> Result<Vec<u8>, CipLabelerError> {
edges
.iter()
.map(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol())
})
.collect()
}
#[test]
fn ciplabeler_sort_priority_orders_edges_and_reports_unique_like_rdkit() {
let priority = CipPriority::new(true, false);
assert!(priority.is_unique());
assert!(!priority.is_pseudo_asymetric());
let molecule = Molecule::from_smiles("C(F)(Cl)Br").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
let rule = AtomicNumberSortRule;
let sorter = CipSort::new(&rule);
let mut context = CipLabelerContext::new(0);
assert_eq!(sorter.get_rules().len(), 1);
let priority = sorter
.prioritize(&mut digraph, &mut context, root, &mut edges, true)
.unwrap();
assert!(priority.is_unique());
assert!(!priority.is_pseudo_asymetric());
assert_eq!(
edge_end_atomic_nums(&digraph, &edges).unwrap(),
vec![35, 17, 9, 1]
);
}
#[test]
fn ciplabeler_sort_priority_groups_equal_edges_like_rdkit() {
let molecule = Molecule::from_smiles("CC(C)N").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
let rule = AtomicNumberSortRule;
let sorter = CipSort::from_rules(vec![&rule]);
let mut context = CipLabelerContext::new(0);
let priority = sorter
.prioritize(&mut digraph, &mut context, root, &mut edges, true)
.unwrap();
let groups = sorter
.get_groups(&mut digraph, &mut context, &edges)
.unwrap();
assert!(!priority.is_unique());
assert!(!priority.is_pseudo_asymetric());
assert_eq!(
edge_end_atomic_nums(&digraph, &edges).unwrap(),
vec![7, 6, 6, 1]
);
assert_eq!(
groups.iter().map(Vec::len).collect::<Vec<_>>(),
vec![1, 2, 1]
);
}
#[test]
fn ciplabeler_sort_priority_counts_single_pseudoasym_comparison_like_rdkit() {
let molecule = Molecule::from_smiles("FON").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
let rule = PseudoAsymSortRule;
let sorter = CipSort::new(&rule);
let mut context = CipLabelerContext::new(0);
let priority = sorter
.prioritize(&mut digraph, &mut context, root, &mut edges, true)
.unwrap();
assert!(priority.is_unique());
assert!(priority.is_pseudo_asymetric());
assert_eq!(edge_end_atomic_nums(&digraph, &edges).unwrap(), vec![9, 7]);
}
#[test]
fn ciplabeler_sequence_rule_get_bond_label_matches_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let rule = BondLabelSortRule;
let real_bond = edges
.iter()
.copied()
.find(|edge_id| digraph.edge(*edge_id).get_bond_idx().is_some())
.unwrap();
let implicit_h = edges
.iter()
.copied()
.find(|edge_id| digraph.edge(*edge_id).get_bond_idx().is_none())
.unwrap();
assert_eq!(
rule.get_bond_label(digraph.edge(real_bond)),
Descriptor::None
);
assert_eq!(
rule.get_bond_label(digraph.edge(implicit_h)),
Descriptor::None
);
digraph.edges[real_bond.index()].set_aux(Descriptor::seqTrans);
digraph.edges[implicit_h.index()].set_aux(Descriptor::seqCis);
assert_eq!(
rule.get_bond_label(digraph.edge(real_bond)),
Descriptor::seqTrans
);
assert_eq!(
rule.get_bond_label(digraph.edge(implicit_h)),
Descriptor::None
);
}
#[test]
fn ciplabeler_sequence_rule_get_comparison_deep_and_shallow_match_rdkit() {
let molecule = Molecule::from_smiles("C(CF)CCl").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let carbon_edges = edges
.iter()
.copied()
.filter(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 6
})
.collect::<Vec<_>>();
assert_eq!(carbon_edges.len(), 2);
let rule = AtomicNumberSortRule;
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.get_comparison(
&mut digraph,
&mut context,
carbon_edges[0],
carbon_edges[1],
false,
)
.unwrap(),
0
);
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.get_comparison(
&mut digraph,
&mut context,
carbon_edges[0],
carbon_edges[1],
true,
)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_sequence_rule_recursive_compare_returns_early_nonzero_like_rdkit() {
let molecule = Molecule::from_smiles("C(F)Cl").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let fluorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 9
})
.unwrap();
let chlorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 17
})
.unwrap();
let rule = AtomicNumberSortRule;
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.recursive_compare(&mut digraph, &mut context, fluorine, chlorine)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_sequence_rule_recursive_compare_enforces_remaining_call_count_like_rdkit() {
let molecule = Molecule::from_smiles("C(F)Cl").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let fluorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 9
})
.unwrap();
let chlorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 17
})
.unwrap();
let rule = AtomicNumberSortRule;
let mut context = CipLabelerContext::with_remaining_call_count(1);
assert!(matches!(
rule.recursive_compare(&mut digraph, &mut context, fluorine, chlorine),
Err(CipLabelerError::MaxIterationsExceeded)
));
}
#[test]
fn ciplabeler_sequence_rule_sort_delegates_to_single_rule_sorter_like_rdkit() {
let molecule = Molecule::from_smiles("C(F)(Cl)Br").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
let rule = AtomicNumberSortRule;
let mut context = CipLabelerContext::new(0);
let priority = rule
.sort(&mut digraph, &mut context, root, &mut edges, true)
.unwrap();
assert!(priority.is_unique());
assert!(!priority.is_pseudo_asymetric());
assert_eq!(
edge_end_atomic_nums(&digraph, &edges).unwrap(),
vec![35, 17, 9, 1]
);
}
#[test]
fn ciplabeler_sequence_rule_are_up_edges_skip_and_error_match_rdkit() {
let molecule = Molecule::from_smiles("CCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(0)
})
.unwrap();
let right = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(2)
})
.unwrap();
let left_node = digraph.edge(left).get_end();
let right_node = digraph.edge(right).get_end();
digraph.node_edges(left_node).unwrap();
let rule = AtomicNumberSortRule;
assert!(
rule.are_up_edges(&digraph, left_node, left_node, left, left)
.unwrap()
);
assert!(matches!(
rule.are_up_edges(&digraph, left_node, root, left, right),
Err(CipLabelerError::UnexpectedUpEdgeOrdering)
));
}
#[test]
fn ciplabeler_rules_constructor_count_and_null_error_match_rdkit() {
let rules = CipRules::new(vec![
Box::new(AlwaysEqualSortRule),
Box::new(AtomicNumberSortRule),
])
.unwrap();
assert_eq!(rules.get_num_sub_rules(), 2);
assert_eq!(rules.get_sorter().get_rules().len(), 1);
let mut empty = CipRules::new(Vec::new()).unwrap();
assert!(matches!(
empty.add(None),
Err(CipLabelerError::NoSequenceRuleProvided)
));
}
#[test]
fn ciplabeler_rules_compare_tries_subrules_in_order_like_rdkit() {
let molecule = Molecule::from_smiles("C(F)Cl").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let fluorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 9
})
.unwrap();
let chlorine = edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 17
})
.unwrap();
let rules = CipRules::new(vec![
Box::new(AlwaysEqualSortRule),
Box::new(AtomicNumberSortRule),
])
.unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
rules
.get_comparison(&mut digraph, &mut context, fluorine, chlorine, false)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_rules_get_comparison_ignores_deep_flag_like_rdkit() {
let molecule = Molecule::from_smiles("C(CF)CCl").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let edges = digraph.node_edges(root).unwrap();
let carbon_edges = edges
.iter()
.copied()
.filter(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).get_atomic_num(digraph.mol()).unwrap() == 6
})
.collect::<Vec<_>>();
assert_eq!(carbon_edges.len(), 2);
let rules = CipRules::new(vec![Box::new(AtomicNumberSortRule)]).unwrap();
let mut shallow_context = CipLabelerContext::new(0);
let mut deep_context = CipLabelerContext::new(0);
let shallow = rules
.get_comparison(
&mut digraph,
&mut shallow_context,
carbon_edges[0],
carbon_edges[1],
false,
)
.unwrap();
let deep = rules
.get_comparison(
&mut digraph,
&mut deep_context,
carbon_edges[0],
carbon_edges[1],
true,
)
.unwrap();
assert_eq!(shallow, -1);
assert_eq!(deep, shallow);
}
#[test]
fn ciplabeler_rules_sort_uses_rules_own_sorter_like_rdkit() {
let molecule = Molecule::from_smiles("C(F)(Cl)Br").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
let rules = CipRules::new(vec![
Box::new(AlwaysEqualSortRule),
Box::new(AtomicNumberSortRule),
])
.unwrap();
let mut context = CipLabelerContext::new(0);
let priority = rules
.sort(&mut digraph, &mut context, root, &mut edges, true)
.unwrap();
assert!(priority.is_unique());
assert!(!priority.is_pseudo_asymetric());
assert_eq!(
edge_end_atomic_nums(&digraph, &edges).unwrap(),
vec![35, 17, 9, 1]
);
}
#[test]
fn ciplabeler_rules_rule1a_compares_fractional_atomic_number_like_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let carbon_node = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
let fractional_node = digraph
.add_node(vec![], Some(1), RationalI32::new(13, 2), 2, 0)
.unwrap();
digraph.add_edge(root, Some(0), carbon_node);
let carbon = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(root, Some(0), fractional_node);
let fractional_n = CipEdgeId::new(digraph.edges.len() - 1);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule1a
.compare(&mut digraph, &mut context, carbon, fractional_n)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_rules_rule1b_orders_ring_duplicates_like_rdkit_non_iupac_2013() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let non_duplicate_node = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
let ring_duplicate_node = digraph
.add_node(
vec![],
Some(1),
RationalI32::new(6, 1),
3,
CipNode::RING_DUPLICATE,
)
.unwrap();
digraph.add_edge(root, Some(0), non_duplicate_node);
let non_duplicate = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(root, Some(0), ring_duplicate_node);
let ring_duplicate = CipEdgeId::new(digraph.edges.len() - 1);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule1b
.compare(&mut digraph, &mut context, ring_duplicate, non_duplicate)
.unwrap(),
1
);
assert_eq!(
CipRule1b
.compare(&mut digraph, &mut context, non_duplicate, ring_duplicate)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_rules_rule2_compares_isotopic_mass_like_rdkit() {
let mut builder = MoleculeBuilder::new();
let center = builder.add_atom(AtomSpec::new(Element::C));
let c12 = builder.add_atom(AtomSpec::new(Element::C).with_isotope(12));
let c13 = builder.add_atom(AtomSpec::new(Element::C).with_isotope(13));
builder
.add_bond(BondSpec::new(center, c12, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(center, c13, BondOrder::Single))
.unwrap();
let molecule = builder.build().unwrap();
let mut digraph = CipDigraph::new(&molecule, center.index(), false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let c13 = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(c13.index())
})
.unwrap();
let c12 = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(c12.index())
})
.unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule2
.compare(&mut digraph, &mut context, c12, c13)
.unwrap(),
-1
);
let mut context = CipLabelerContext::new(0);
let c13_mass = digraph.node(digraph.edge(c13).get_end()).get_atomic_mass();
let c12_mass = digraph.node(digraph.edge(c12).get_end()).get_atomic_mass();
assert!(c13_mass > c12_mass);
assert_eq!(
CipRule2
.compare(&mut digraph, &mut context, c13, c12)
.unwrap(),
1
);
}
#[test]
fn ciplabeler_rules_rule3_orders_e_z_aux_labels_like_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let mut edges = digraph.node_edges(root).unwrap();
assert!(edges.len() >= 2);
let e_edge = edges[0];
let z_edge = edges[1];
let e_node = digraph.edge(e_edge).get_end();
let z_node = digraph.edge(z_edge).get_end();
digraph.nodes[e_node.index()].set_aux(Descriptor::E);
digraph.nodes[z_node.index()].set_aux(Descriptor::Z);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule3
.compare(&mut digraph, &mut context, e_edge, z_edge)
.unwrap(),
-1
);
edges.swap(0, 1);
assert_eq!(
CipRule3
.compare(&mut digraph, &mut context, edges[0], edges[1])
.unwrap(),
1
);
}
#[test]
fn ciplabeler_pairlist_ref_and_add_match_rdkit() {
assert_eq!(CipPairList::ref_descriptor(Descriptor::R), Descriptor::R);
assert_eq!(CipPairList::ref_descriptor(Descriptor::M), Descriptor::R);
assert_eq!(
CipPairList::ref_descriptor(Descriptor::seqCis),
Descriptor::R
);
assert_eq!(CipPairList::ref_descriptor(Descriptor::S), Descriptor::S);
assert_eq!(CipPairList::ref_descriptor(Descriptor::P), Descriptor::S);
assert_eq!(
CipPairList::ref_descriptor(Descriptor::seqTrans),
Descriptor::S
);
assert_eq!(
CipPairList::ref_descriptor(Descriptor::None),
Descriptor::None
);
assert_eq!(CipPairList::ref_descriptor(Descriptor::r), Descriptor::None);
let mut pairs = CipPairList::new();
assert!(!pairs.add(Descriptor::None));
assert!(!pairs.add(Descriptor::r));
assert!(pairs.add(Descriptor::M));
assert_eq!(pairs.get_ref_descriptor(), Descriptor::R);
assert_eq!(pairs.to_rdkit_string(), "R:");
}
#[test]
fn ciplabeler_pairlist_pairing_bits_and_string_match_rdkit() {
let mut pairs = CipPairList::with_ref(Descriptor::R);
pairs.add(Descriptor::R);
pairs.add(Descriptor::P);
pairs.add(Descriptor::R);
assert_eq!(pairs.get_pairing(), (1_u64 << 62) | (1_u64 << 60));
assert_eq!(pairs.to_rdkit_string(), "R:lul");
}
#[test]
fn ciplabeler_pairlist_head_tail_constructor_filters_descriptors_like_rdkit() {
let mut head = CipPairList::with_ref(Descriptor::R);
head.add(Descriptor::S);
let mut tail = CipPairList::new();
tail.add(Descriptor::None);
tail.add(Descriptor::seqCis);
tail.add(Descriptor::Unknown);
let combined = CipPairList::from_head_tail(&head, &tail);
assert_eq!(combined.to_rdkit_string(), "R:ul");
assert_eq!(combined.get_pairing(), 1_u64 << 61);
}
#[test]
fn ciplabeler_pairlist_compare_to_and_order_match_rdkit() {
let mut like = CipPairList::with_ref(Descriptor::R);
like.add(Descriptor::R);
let mut unlike = CipPairList::with_ref(Descriptor::R);
unlike.add(Descriptor::S);
assert_eq!(like.compare_to(&unlike).unwrap(), 1);
assert_eq!(unlike.compare_to(&like).unwrap(), -1);
assert!(unlike < like);
let mut shorter = CipPairList::with_ref(Descriptor::R);
shorter.add(Descriptor::R);
let mut longer = shorter.clone();
longer.add(Descriptor::R);
assert!(matches!(
shorter.compare_to(&longer),
Err(CipLabelerError::DescriptorListLengthMismatch)
));
}
#[test]
fn ciplabeler_rules_rule4a_orders_descriptor_classes_like_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let node_r = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
let node_s = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
digraph.add_edge(root, Some(0), node_r);
let edge_r = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(root, Some(0), node_s);
let edge_s = CipEdgeId::new(digraph.edges.len() - 1);
digraph.edges[edge_r.index()].set_aux(Descriptor::R);
digraph.edges[edge_s.index()].set_aux(Descriptor::r);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule4a
.compare(&mut digraph, &mut context, edge_r, edge_s)
.unwrap(),
1
);
digraph.edges[edge_r.index()].set_aux(Descriptor::None);
digraph.edges[edge_s.index()].set_aux(Descriptor::None);
digraph.nodes[node_r.index()].set_aux(Descriptor::R);
digraph.nodes[node_s.index()].set_aux(Descriptor::r);
assert_eq!(
CipRule4a
.compare(&mut digraph, &mut context, edge_r, edge_s)
.unwrap(),
1
);
digraph.nodes[node_r.index()].set_aux(Descriptor::SP_4);
assert!(matches!(
CipRule4a.compare(&mut digraph, &mut context, edge_r, edge_s),
Err(CipLabelerError::InvalidStereoDescriptor)
));
}
#[test]
fn ciplabeler_rules_rule4c_orders_lowercase_descriptor_classes_like_rdkit() {
let molecule = Molecule::from_smiles("CC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 0, false).unwrap();
let root = digraph.get_current_root();
let node_m = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
let node_p = digraph
.add_node(vec![], Some(1), RationalI32::new(6, 1), 2, 0)
.unwrap();
digraph.add_edge(root, Some(0), node_m);
let edge_m = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(root, Some(0), node_p);
let edge_p = CipEdgeId::new(digraph.edges.len() - 1);
digraph.edges[edge_m.index()].set_aux(Descriptor::m);
digraph.edges[edge_p.index()].set_aux(Descriptor::p);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule4c
.compare(&mut digraph, &mut context, edge_m, edge_p)
.unwrap(),
1
);
digraph.edges[edge_m.index()].set_aux(Descriptor::None);
digraph.edges[edge_p.index()].set_aux(Descriptor::None);
digraph.nodes[node_m.index()].set_aux(Descriptor::r);
digraph.nodes[node_p.index()].set_aux(Descriptor::s);
assert_eq!(
CipRule4c
.compare(&mut digraph, &mut context, edge_m, edge_p)
.unwrap(),
1
);
}
#[test]
fn ciplabeler_rules_rule4b_nonroot_ref_branch_matches_rdkit() {
let molecule = Molecule::from_smiles("CCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let branch = digraph
.add_node(vec![1, 1, 0], Some(0), RationalI32::new(6, 1), 2, 0)
.unwrap();
digraph.add_edge(root, Some(0), branch);
let r_child = digraph
.add_node(vec![], Some(0), RationalI32::new(6, 1), 3, 0)
.unwrap();
let s_child = digraph
.add_node(vec![], Some(0), RationalI32::new(6, 1), 3, 0)
.unwrap();
digraph.nodes[r_child.index()].set_aux(Descriptor::R);
digraph.nodes[s_child.index()].set_aux(Descriptor::S);
digraph.add_edge(branch, Some(0), r_child);
let r_edge = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(branch, Some(0), s_child);
let s_edge = CipEdgeId::new(digraph.edges.len() - 1);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule4b::with_ref(Descriptor::R)
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
1
);
assert_eq!(
CipRule4b::with_ref(Descriptor::S)
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
-1
);
assert_eq!(
CipRule4b::new()
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
0
);
}
#[test]
fn ciplabeler_rules_rule4b_reference_descriptor_search_matches_rdkit() {
let molecule = Molecule::from_smiles("CCCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(0)
})
.unwrap();
let right = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(2)
})
.unwrap();
let left_node = digraph.edge(left).get_end();
let right_node = digraph.edge(right).get_end();
digraph.nodes[left_node.index()].set_aux(Descriptor::R);
digraph.nodes[right_node.index()].set_aux(Descriptor::S);
let rule = CipRule4b::new();
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.get_reference_descriptors(None, &mut digraph, &mut context, left_node)
.unwrap(),
vec![Descriptor::R]
);
assert_eq!(
rule.get_reference_descriptors(None, &mut digraph, &mut context, right_node)
.unwrap(),
vec![Descriptor::S]
);
assert!(rule.has_descriptors(&mut digraph, left_node).unwrap());
}
#[test]
fn ciplabeler_rules_rule4b_root_pair_comparison_matches_rdkit() {
let molecule = Molecule::from_smiles("CCCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(0)
})
.unwrap();
let right = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(2)
})
.unwrap();
let left_node = digraph.edge(left).get_end();
let right_node = digraph.edge(right).get_end();
digraph.nodes[left_node.index()].set_aux(Descriptor::R);
digraph.nodes[right_node.index()].set_aux(Descriptor::S);
let rule = CipRule4b::new();
let sort_rules = vec![&rule as &dyn CipSequenceRule];
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.compare_with_sort_rules(
Some(&sort_rules),
&mut digraph,
&mut context,
left,
right,
)
.unwrap(),
0
);
}
#[test]
fn ciplabeler_rules_rule5new_rule6_nonroot_ref_branch_matches_rdkit() {
let molecule = Molecule::from_smiles("CCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let branch = digraph
.add_node(vec![1, 1, 0], Some(0), RationalI32::new(6, 1), 2, 0)
.unwrap();
digraph.add_edge(root, Some(0), branch);
let r_child = digraph
.add_node(vec![], Some(0), RationalI32::new(6, 1), 3, 0)
.unwrap();
let s_child = digraph
.add_node(vec![], Some(0), RationalI32::new(6, 1), 3, 0)
.unwrap();
digraph.nodes[r_child.index()].set_aux(Descriptor::R);
digraph.nodes[s_child.index()].set_aux(Descriptor::S);
digraph.add_edge(branch, Some(0), r_child);
let r_edge = CipEdgeId::new(digraph.edges.len() - 1);
digraph.add_edge(branch, Some(0), s_child);
let s_edge = CipEdgeId::new(digraph.edges.len() - 1);
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule5New::with_ref(Descriptor::R)
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
1
);
assert_eq!(
CipRule5New::with_ref(Descriptor::S)
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
-1
);
assert_eq!(
CipRule5New::new()
.compare(&mut digraph, &mut context, r_edge, s_edge)
.unwrap(),
0
);
}
#[test]
fn ciplabeler_rules_rule5new_rule6_root_pair_comparison_matches_rdkit() {
let molecule = Molecule::from_smiles("CCCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(0)
})
.unwrap();
let right = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(2)
})
.unwrap();
let left_node = digraph.edge(left).get_end();
let right_node = digraph.edge(right).get_end();
digraph.nodes[left_node.index()].set_aux(Descriptor::R);
digraph.nodes[right_node.index()].set_aux(Descriptor::S);
let rule = CipRule5New::new();
let sort_rules = vec![&rule as &dyn CipSequenceRule];
let mut context = CipLabelerContext::new(0);
assert_eq!(
rule.compare_with_sort_rules(
Some(&sort_rules),
&mut digraph,
&mut context,
left,
right,
)
.unwrap(),
2
);
}
#[test]
fn ciplabeler_rules_rule5new_rule6_rule6_ref_atom_matches_rdkit() {
let molecule = Molecule::from_smiles("CCC").unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_current_root();
let root_edges = digraph.node_edges(root).unwrap();
let left = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(0)
})
.unwrap();
let right = root_edges
.iter()
.copied()
.find(|edge_id| {
let end = digraph.edge(*edge_id).get_end();
digraph.node(end).atom_idx() == Some(2)
})
.unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
CipRule6
.compare(&mut digraph, &mut context, left, right)
.unwrap(),
0
);
digraph.set_rule6_ref(Some(0)).unwrap();
assert_eq!(
CipRule6
.compare(&mut digraph, &mut context, left, right)
.unwrap(),
1
);
digraph.set_rule6_ref(Some(2)).unwrap();
assert_eq!(
CipRule6
.compare(&mut digraph, &mut context, left, right)
.unwrap(),
-1
);
}
#[test]
fn ciplabeler_configuration_constructors_accessors_and_carriers_match_rdkit() {
let molecule = Molecule::from_smiles("CCO").unwrap();
let mut config = CipConfiguration::new(&molecule, 1).unwrap();
assert_eq!(CipConfiguration::IMPLICIT_H, usize::MAX);
assert_eq!(config.get_focus(), 1);
assert_eq!(config.get_foci(), &[1]);
assert!(config.get_carriers().is_empty());
assert_eq!(config.get_digraph().get_original_root(), CipNodeId::new(0));
assert_eq!(
config.get_digraph().node(CipNodeId::new(0)).atom_idx(),
Some(1)
);
config.set_carriers(vec![Some(0), Some(2), None]);
assert_eq!(config.get_carriers(), &[Some(0), Some(2), None]);
let mut multi = CipConfiguration::with_foci(&molecule, vec![0, 1], true).unwrap();
assert_eq!(multi.get_focus(), 0);
assert_eq!(multi.get_foci(), &[0, 1]);
assert_eq!(multi.get_digraph().get_original_root(), CipNodeId::new(0));
assert!(matches!(
CipConfiguration::with_foci(&molecule, Vec::new(), false),
Err(CipLabelerError::EmptyConfigurationFoci)
));
}
#[test]
fn ciplabeler_configuration_parity4_matches_rdkit_table() {
let reference = [0, 1, 2, 3];
let even = [
[0, 1, 2, 3],
[0, 2, 3, 1],
[0, 3, 1, 2],
[1, 0, 3, 2],
[1, 2, 0, 3],
[1, 3, 2, 0],
[2, 0, 1, 3],
[2, 1, 3, 0],
[2, 3, 0, 1],
[3, 0, 2, 1],
[3, 1, 0, 2],
[3, 2, 1, 0],
];
let odd = [
[0, 1, 3, 2],
[0, 2, 1, 3],
[0, 3, 2, 1],
[1, 0, 2, 3],
[1, 2, 3, 0],
[1, 3, 0, 2],
[2, 0, 3, 1],
[2, 1, 0, 3],
[2, 3, 1, 0],
[3, 0, 1, 2],
[3, 1, 2, 0],
[3, 2, 0, 1],
];
for target in even {
assert_eq!(
CipConfiguration::parity4(&target, &reference).unwrap(),
2,
"target={target:?}"
);
}
for target in odd {
assert_eq!(
CipConfiguration::parity4(&target, &reference).unwrap(),
1,
"target={target:?}"
);
}
assert_eq!(
CipConfiguration::parity4(&[9, 8, 7, 6], &reference).unwrap(),
0
);
assert!(matches!(
CipConfiguration::parity4(&[0, 1, 2], &reference),
Err(CipLabelerError::ParityVectorsMustHaveSize4)
));
}
#[test]
fn ciplabeler_configuration_base_label_returns_unknown_like_rdkit() {
let molecule = Molecule::from_smiles("CCO").unwrap();
let config = CipConfiguration::new(&molecule, 1).unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let node = digraph.get_current_root();
let rules = CipRules::new(vec![Box::new(CipRule1a)]).unwrap();
assert_eq!(
config.label(node, &mut digraph, &rules),
Descriptor::Unknown
);
}
fn tetrahedral_molecule(tag: ChiralTag) -> Molecule {
let mut builder = MoleculeBuilder::new();
let center = builder.add_atom(AtomSpec::new(Element::C).with_chiral_tag(tag));
let fluorine = builder.add_atom(AtomSpec::new(Element::F));
let chlorine = builder.add_atom(AtomSpec::new(Element::CL));
let bromine = builder.add_atom(AtomSpec::new(Element::BR));
let iodine = builder.add_atom(AtomSpec::new(Element::I));
builder
.add_bond(BondSpec::new(center, fluorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(center, chlorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(center, bromine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(center, iodine, BondOrder::Single))
.unwrap();
builder.build().unwrap()
}
#[test]
fn ciplabeler_tetrahedral_constructor_builds_carriers_like_rdkit() {
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let config = CipTetrahedral::new(&molecule, 0).unwrap();
assert_eq!(config.get_focus(), 0);
assert_eq!(config.get_carriers(), &[Some(1), Some(2), Some(3), Some(4)]);
let implicit_h = Molecule::from_smiles_with_sanitize("C[C@H](F)Cl", false).unwrap();
let implicit_config = CipTetrahedral::new(&implicit_h, 1).unwrap();
assert_eq!(
implicit_config.get_carriers(),
&[Some(0), Some(2), Some(3), Some(1)]
);
let methane = Molecule::from_smiles("C").unwrap();
assert!(matches!(
CipTetrahedral::new(&methane, 0),
Err(CipLabelerError::BadTetrahedralConfig)
));
}
#[test]
fn ciplabeler_tetrahedral_primary_label_records_atom_props_like_rdkit() {
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let mut config = CipTetrahedral::new(&molecule, 0).unwrap();
assert!(!config.has_primary_label());
config.ranked_anchors = vec![4, 3, 2, 1];
config.set_primary_label(Descriptor::S).unwrap();
assert!(config.has_primary_label());
assert_eq!(
config.primary_label(),
Some(&super::CipAtomPrimaryLabel {
atom_idx: 0,
cip_code: "S",
cip_neighbor_order: vec![4, 3, 2, 1],
})
);
config.reset_primary_label();
assert!(!config.has_primary_label());
assert!(matches!(
config.set_primary_label(Descriptor::E),
Err(CipLabelerError::DescriptorNotSupportedForAtoms)
));
assert!(matches!(
config.set_primary_label(Descriptor::Unknown),
Err(CipLabelerError::InvalidAtomDescriptor)
));
}
#[test]
fn ciplabeler_tetrahedral_label_assigns_r_s_and_ranked_anchors_like_rdkit() {
let rules = cip_all_rules().unwrap();
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let mut config = CipTetrahedral::new(&molecule, 0).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::S);
assert_eq!(config.ranked_anchors(), &[4, 3, 2, 1]);
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCw);
let mut config = CipTetrahedral::new(&molecule, 0).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::R);
assert_eq!(config.ranked_anchors(), &[4, 3, 2, 1]);
}
#[test]
fn ciplabeler_tetrahedral_label_external_digraph_changes_root_like_rdkit() {
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let mut config = CipTetrahedral::new(&molecule, 0).unwrap();
let mut external = CipDigraph::new(&molecule, 0, false).unwrap();
let root = external.get_original_root();
let rules = cip_all_rules().unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config
.label_with_external_digraph(root, &mut external, &rules, &mut context)
.unwrap(),
Descriptor::S
);
assert_eq!(external.get_current_root(), root);
}
#[test]
fn ciplabeler_tetrahedral_label_returns_ns_or_unknown_for_unresolved_cases_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("C[C@H](C)C", false).unwrap();
let mut config = CipTetrahedral::new(&molecule, 1).unwrap();
let rules = cip_all_rules().unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config.label(&rules, &mut context).unwrap(),
Descriptor::Unknown
);
let leaf_molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let mut config = CipTetrahedral::new(&leaf_molecule, 0).unwrap();
let mut external = CipDigraph::new(&leaf_molecule, 0, false).unwrap();
let root = external.get_original_root();
let leaf_edge = external.node_edges(root).unwrap()[0];
let leaf = external.edge(leaf_edge).get_end();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config
.label_with_external_digraph(leaf, &mut external, &rules, &mut context)
.unwrap(),
Descriptor::ns
);
}
fn sp2_bond_molecule(stereo: BondStereo, stereo_atoms: (usize, usize)) -> Molecule {
let mut builder = MoleculeBuilder::new();
let fluorine = builder.add_atom(AtomSpec::new(Element::F));
let begin = builder.add_atom(AtomSpec::new(Element::C));
let end = builder.add_atom(AtomSpec::new(Element::C));
let chlorine = builder.add_atom(AtomSpec::new(Element::CL));
let bromine = builder.add_atom(AtomSpec::new(Element::BR));
let iodine = builder.add_atom(AtomSpec::new(Element::I));
builder
.add_bond(
BondSpec::new(begin, end, BondOrder::Double)
.with_stereo(stereo)
.with_stereo_atoms(
match stereo_atoms.0 {
0 => fluorine,
3 => chlorine,
4 => bromine,
5 => iodine,
_ => panic!("unexpected stereo atom fixture index"),
},
match stereo_atoms.1 {
0 => fluorine,
3 => chlorine,
4 => bromine,
5 => iodine,
_ => panic!("unexpected stereo atom fixture index"),
},
),
)
.unwrap();
builder
.add_bond(BondSpec::new(begin, fluorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(begin, bromine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, chlorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, iodine, BondOrder::Single))
.unwrap();
builder.build().unwrap()
}
fn sp2_ez_bond_molecule_without_explicit_stereo_atoms() -> Molecule {
let mut builder = MoleculeBuilder::new();
let fluorine = builder.add_atom(AtomSpec::new(Element::F).with_prop("_CIPRank", "10"));
let begin = builder.add_atom(AtomSpec::new(Element::C));
let end = builder.add_atom(AtomSpec::new(Element::C));
let chlorine = builder.add_atom(AtomSpec::new(Element::CL).with_prop("_CIPRank", "20"));
let bromine = builder.add_atom(AtomSpec::new(Element::BR).with_prop("_CIPRank", "30"));
let iodine = builder.add_atom(AtomSpec::new(Element::I).with_prop("_CIPRank", "40"));
builder
.add_bond(BondSpec::new(begin, end, BondOrder::Double).with_stereo(BondStereo::E))
.unwrap();
builder
.add_bond(BondSpec::new(begin, fluorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(begin, bromine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, chlorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, iodine, BondOrder::Single))
.unwrap();
builder.build().unwrap()
}
#[test]
fn ciplabeler_sp2bond_constructor_builds_carriers_like_rdkit_find_stereo_atoms() {
let molecule = sp2_bond_molecule(BondStereo::Cis, (4, 5));
let config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Cis).unwrap();
assert_eq!(config.get_foci(), &[1, 2]);
assert_eq!(config.get_carriers(), &[Some(4), Some(5)]);
assert!(matches!(
CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::E),
Err(CipLabelerError::BadSp2BondConfig)
));
assert!(matches!(
CipSp2Bond::new(&molecule, 1, 1, 0, BondStereo::Cis),
Err(CipLabelerError::BadSp2BondFoci)
));
let mut builder = MoleculeBuilder::new();
let a0 = builder.add_atom(AtomSpec::new(Element::C));
let a1 = builder.add_atom(AtomSpec::new(Element::C));
builder
.add_bond(BondSpec::new(a0, a1, BondOrder::Double).with_stereo(BondStereo::E))
.unwrap();
let missing = builder.build().unwrap();
assert!(matches!(
CipSp2Bond::new(&missing, 0, 0, 1, BondStereo::Cis),
Err(CipLabelerError::IncorrectNumberOfStereoAtoms)
));
}
#[test]
fn ciplabeler_sp2bond_constructor_finds_highest_cip_neighbors_like_rdkit() {
let molecule = sp2_ez_bond_molecule_without_explicit_stereo_atoms();
let config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Trans).unwrap();
assert_eq!(config.get_foci(), &[1, 2]);
assert_eq!(config.get_carriers(), &[Some(4), Some(5)]);
}
#[test]
fn ciplabeler_sp2bond_primary_label_records_bond_props_like_rdkit() {
let molecule = sp2_bond_molecule(BondStereo::Cis, (4, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Cis).unwrap();
assert!(!config.has_primary_label());
config.ranked_anchors = vec![4, 5];
config.set_primary_label(Descriptor::Z).unwrap();
assert!(config.has_primary_label());
assert_eq!(
config.primary_label(),
Some(&super::CipBondPrimaryLabel {
bond_idx: 0,
stereo_atoms: [4, 5],
stereo: BondStereo::Cis,
cip_code: "Z",
cip_neighbor_order: vec![4, 5],
})
);
config.reset_primary_label();
assert!(!config.has_primary_label());
assert!(matches!(
config.set_primary_label(Descriptor::R),
Err(CipLabelerError::DescriptorNotSupportedForDoubleBonds)
));
assert!(matches!(
config.set_primary_label(Descriptor::Unknown),
Err(CipLabelerError::InvalidBondDescriptor)
));
}
#[test]
fn ciplabeler_sp2bond_label_assigns_e_z_and_ranked_anchors_like_rdkit() {
let rules = cip_all_rules().unwrap();
let molecule = sp2_bond_molecule(BondStereo::Cis, (4, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Cis).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::Z);
assert_eq!(config.ranked_anchors(), &[4, 5]);
let molecule = sp2_bond_molecule(BondStereo::Trans, (4, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Trans).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::E);
assert_eq!(config.ranked_anchors(), &[4, 5]);
}
#[test]
fn ciplabeler_sp2bond_label_flips_config_when_carrier_is_not_top_priority_like_rdkit() {
let rules = cip_all_rules().unwrap();
let molecule = sp2_bond_molecule(BondStereo::Cis, (0, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Cis).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::E);
assert_eq!(config.ranked_anchors(), &[4, 5]);
let molecule = sp2_bond_molecule(BondStereo::Trans, (0, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Trans).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::Z);
assert_eq!(config.ranked_anchors(), &[4, 5]);
}
#[test]
fn ciplabeler_sp2bond_label_external_digraph_matches_rdkit_overload() {
let molecule = sp2_bond_molecule(BondStereo::Cis, (4, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, BondStereo::Cis).unwrap();
let mut external = CipDigraph::new(&molecule, 1, false).unwrap();
let root = external.get_original_root();
let rules = cip_all_rules().unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config
.label_with_external_digraph(root, &mut external, &rules, &mut context)
.unwrap(),
Descriptor::Z
);
assert_eq!(config.ranked_anchors(), &[4, 5]);
}
fn atropisomer_bond_molecule(stereo: BondStereo, low_index_carriers: bool) -> Molecule {
atropisomer_bond_molecule_with_carrier_positions(
stereo,
low_index_carriers,
low_index_carriers,
)
}
fn atropisomer_bond_molecule_with_carrier_positions(
stereo: BondStereo,
begin_low_index_carrier: bool,
end_low_index_carrier: bool,
) -> Molecule {
let mut builder = MoleculeBuilder::new();
let fluorine = builder.add_atom(AtomSpec::new(Element::F));
let begin = builder.add_atom(AtomSpec::new(Element::C));
let end = builder.add_atom(AtomSpec::new(Element::C));
let chlorine = builder.add_atom(AtomSpec::new(Element::CL));
let bromine = builder.add_atom(AtomSpec::new(Element::BR));
let iodine = builder.add_atom(AtomSpec::new(Element::I));
builder
.add_bond(BondSpec::new(begin, end, BondOrder::Single).with_stereo(stereo))
.unwrap();
if begin_low_index_carrier {
builder
.add_bond(BondSpec::new(begin, fluorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(begin, bromine, BondOrder::Single))
.unwrap();
} else {
builder
.add_bond(BondSpec::new(begin, bromine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(begin, fluorine, BondOrder::Single))
.unwrap();
}
if end_low_index_carrier {
builder
.add_bond(BondSpec::new(end, chlorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, iodine, BondOrder::Single))
.unwrap();
} else {
builder
.add_bond(BondSpec::new(end, iodine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, chlorine, BondOrder::Single))
.unwrap();
}
builder.build().unwrap()
}
fn atropisomer_bond_molecule_one_side_second_carrier(stereo: BondStereo) -> Molecule {
let mut builder = MoleculeBuilder::new();
let bromine = builder.add_atom(AtomSpec::new(Element::BR));
let begin = builder.add_atom(AtomSpec::new(Element::C));
let end = builder.add_atom(AtomSpec::new(Element::C));
let chlorine = builder.add_atom(AtomSpec::new(Element::CL));
let fluorine = builder.add_atom(AtomSpec::new(Element::F));
let iodine = builder.add_atom(AtomSpec::new(Element::I));
builder
.add_bond(BondSpec::new(begin, end, BondOrder::Single).with_stereo(stereo))
.unwrap();
builder
.add_bond(BondSpec::new(begin, bromine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(begin, fluorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, chlorine, BondOrder::Single))
.unwrap();
builder
.add_bond(BondSpec::new(end, iodine, BondOrder::Single))
.unwrap();
builder.build().unwrap()
}
#[test]
fn ciplabeler_atropisomerbond_constructor_builds_carriers_like_rdkit() {
let molecule = atropisomer_bond_molecule(BondStereo::AtropCcw, true);
let config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCcw).unwrap();
assert_eq!(config.get_foci(), &[1, 2]);
assert_eq!(config.get_carriers(), &[Some(0), Some(3)]);
assert!(matches!(
CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::Cis),
Err(CipLabelerError::BadAtropisomerBondConfig)
));
assert!(matches!(
CipAtropisomerBond::new(&molecule, 1, 1, 2, BondStereo::AtropCcw),
Err(CipLabelerError::BadAtropisomerBondFoci)
));
let mut builder = MoleculeBuilder::new();
let a0 = builder.add_atom(AtomSpec::new(Element::C));
let a1 = builder.add_atom(AtomSpec::new(Element::C));
builder
.add_bond(BondSpec::new(a0, a1, BondOrder::Single).with_stereo(BondStereo::AtropCcw))
.unwrap();
let missing = builder.build().unwrap();
let config = CipAtropisomerBond::new(&missing, 0, 0, 1, BondStereo::AtropCcw).unwrap();
assert!(config.get_carriers().is_empty());
}
#[test]
fn ciplabeler_atropisomerbond_primary_label_records_bond_props_like_rdkit() {
let molecule = atropisomer_bond_molecule(BondStereo::AtropCcw, true);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCcw).unwrap();
assert!(!config.has_primary_label());
config.ranked_anchors = vec![4, 5];
config.set_primary_label(Descriptor::M).unwrap();
assert!(config.has_primary_label());
assert_eq!(
config.primary_label(),
Some(&super::CipAtropisomerBondPrimaryLabel {
bond_idx: 0,
cip_code: "M",
cip_neighbor_order: vec![4, 5],
})
);
config.reset_primary_label();
assert!(!config.has_primary_label());
assert!(matches!(
config.set_primary_label(Descriptor::Z),
Err(CipLabelerError::DescriptorNotSupportedForAtropisomerBonds)
));
assert!(matches!(
config.set_primary_label(Descriptor::Unknown),
Err(CipLabelerError::InvalidBondDescriptor)
));
}
#[test]
fn ciplabeler_atropisomerbond_label_assigns_m_p_and_ranked_anchors_like_rdkit() {
let rules = cip_all_rules().unwrap();
let molecule = atropisomer_bond_molecule(BondStereo::AtropCcw, false);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCcw).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::M);
assert_eq!(config.ranked_anchors(), &[4, 5]);
let molecule = atropisomer_bond_molecule(BondStereo::AtropCw, false);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCw).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::P);
assert_eq!(config.ranked_anchors(), &[4, 5]);
}
#[test]
fn ciplabeler_atropisomerbond_label_flips_config_when_carrier_is_second_like_rdkit() {
let rules = cip_all_rules().unwrap();
let molecule = atropisomer_bond_molecule_one_side_second_carrier(BondStereo::AtropCcw);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCcw).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::P);
assert_eq!(config.ranked_anchors(), &[0, 5]);
let molecule = atropisomer_bond_molecule_one_side_second_carrier(BondStereo::AtropCw);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCw).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(config.label(&rules, &mut context).unwrap(), Descriptor::M);
assert_eq!(config.ranked_anchors(), &[0, 5]);
}
#[test]
fn ciplabeler_atropisomerbond_label_external_digraph_matches_rdkit_overload() {
let molecule = atropisomer_bond_molecule(BondStereo::AtropCcw, false);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, BondStereo::AtropCcw).unwrap();
let mut external = CipDigraph::new(&molecule, 1, true).unwrap();
let root = external.get_original_root();
let rules = cip_all_rules().unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config
.label_with_external_digraph(root, &mut external, &rules, &mut context)
.unwrap(),
Descriptor::M
);
assert_eq!(config.ranked_anchors(), &[4, 5]);
}
#[test]
fn ciplabeler_configs_golden_configuration_parity_and_base_label_match_rdkit() {
let reference = [0, 1, 2, 3];
let rdkit_parity4_golden = [
([0, 1, 2, 3], 2),
([0, 1, 3, 2], 1),
([0, 2, 1, 3], 1),
([0, 2, 3, 1], 2),
([0, 3, 1, 2], 2),
([0, 3, 2, 1], 1),
([1, 0, 2, 3], 1),
([1, 0, 3, 2], 2),
([1, 2, 0, 3], 2),
([1, 2, 3, 0], 1),
([1, 3, 0, 2], 1),
([1, 3, 2, 0], 2),
([2, 0, 1, 3], 2),
([2, 0, 3, 1], 1),
([2, 1, 0, 3], 1),
([2, 1, 3, 0], 2),
([2, 3, 0, 1], 2),
([2, 3, 1, 0], 1),
([3, 0, 1, 2], 1),
([3, 0, 2, 1], 2),
([3, 1, 0, 2], 2),
([3, 1, 2, 0], 1),
([3, 2, 0, 1], 1),
([3, 2, 1, 0], 2),
];
for (target, expected) in rdkit_parity4_golden {
assert_eq!(
CipConfiguration::parity4(&target, &reference).unwrap(),
expected,
"RDKit parity4 golden mismatch for target {target:?}"
);
}
let molecule = Molecule::from_smiles("CCO").unwrap();
let config = CipConfiguration::new(&molecule, 1).unwrap();
let mut digraph = CipDigraph::new(&molecule, 1, false).unwrap();
let root = digraph.get_original_root();
let rules = cip_all_rules().unwrap();
assert_eq!(
config.label(root, &mut digraph, &rules),
Descriptor::Unknown
);
}
#[test]
fn ciplabeler_configs_golden_tetrahedral_labels_reset_and_neighbor_order_match_rdkit() {
let rules = cip_all_rules().unwrap();
let pseudo_rules = CipRules::new(vec![Box::new(SinglePseudoAsymSortRule)]).unwrap();
let tetrahedral_cases = [
(ChiralTag::TetrahedralCcw, &rules, Descriptor::S, "S"),
(ChiralTag::TetrahedralCw, &rules, Descriptor::R, "R"),
(ChiralTag::TetrahedralCcw, &pseudo_rules, Descriptor::s, "s"),
(ChiralTag::TetrahedralCw, &pseudo_rules, Descriptor::r, "r"),
];
for (tag, rules, expected_descriptor, expected_code) in tetrahedral_cases {
let molecule = tetrahedral_molecule(tag);
let mut config = CipTetrahedral::new(&molecule, 0).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config.label(rules, &mut context).unwrap(),
expected_descriptor
);
assert_eq!(config.ranked_anchors(), &[4, 3, 2, 1]);
config.set_primary_label(expected_descriptor).unwrap();
assert_eq!(
config.primary_label(),
Some(&super::CipAtomPrimaryLabel {
atom_idx: 0,
cip_code: expected_code,
cip_neighbor_order: vec![4, 3, 2, 1],
})
);
assert!(config.has_primary_label());
config.reset_primary_label();
assert_eq!(config.primary_label(), None);
assert!(!config.has_primary_label());
}
}
#[test]
fn ciplabeler_configs_golden_sp2bond_labels_reset_and_neighbor_order_match_rdkit() {
let rules = cip_all_rules().unwrap();
let sp2_cases = [
(BondStereo::Cis, Descriptor::Z, "Z"),
(BondStereo::Trans, Descriptor::E, "E"),
];
for (stereo, expected_descriptor, expected_code) in sp2_cases {
let molecule = sp2_bond_molecule(stereo, (4, 5));
let mut config = CipSp2Bond::new(&molecule, 0, 1, 2, stereo).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config.label(&rules, &mut context).unwrap(),
expected_descriptor
);
assert_eq!(config.ranked_anchors(), &[4, 5]);
config.set_primary_label(expected_descriptor).unwrap();
assert_eq!(
config.primary_label(),
Some(&super::CipBondPrimaryLabel {
bond_idx: 0,
stereo_atoms: [4, 5],
stereo,
cip_code: expected_code,
cip_neighbor_order: vec![4, 5],
})
);
assert!(config.has_primary_label());
config.reset_primary_label();
assert_eq!(config.primary_label(), None);
assert!(!config.has_primary_label());
}
}
#[test]
fn ciplabeler_configs_golden_atropisomerbond_labels_reset_and_neighbor_order_match_rdkit() {
let rules = cip_all_rules().unwrap();
let atrop_cases = [
(BondStereo::AtropCcw, Descriptor::M, "M"),
(BondStereo::AtropCw, Descriptor::P, "P"),
];
for (stereo, expected_descriptor, expected_code) in atrop_cases {
let molecule = atropisomer_bond_molecule(stereo, false);
let mut config = CipAtropisomerBond::new(&molecule, 0, 1, 2, stereo).unwrap();
let mut context = CipLabelerContext::new(0);
assert_eq!(
config.label(&rules, &mut context).unwrap(),
expected_descriptor
);
assert_eq!(config.ranked_anchors(), &[4, 5]);
config.set_primary_label(expected_descriptor).unwrap();
assert_eq!(
config.primary_label(),
Some(&super::CipAtropisomerBondPrimaryLabel {
bond_idx: 0,
cip_code: expected_code,
cip_neighbor_order: vec![4, 5],
})
);
assert!(config.has_primary_label());
config.reset_primary_label();
assert_eq!(config.primary_label(), None);
assert!(!config.has_primary_label());
}
}
#[test]
fn ciplabeler_assign_writes_tetrahedral_labels_and_molecule_computed_prop_like_rdkit() {
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let labeled = assign_cip_labels(&molecule, 0).unwrap();
assert_eq!(molecule.prop("_CIPComputed"), None);
assert_eq!(molecule.atoms()[0].prop("_CIPCode"), None);
assert_eq!(labeled.prop("_CIPComputed"), Some("1"));
assert_eq!(labeled.atoms()[0].prop("_CIPCode"), Some("S"));
assert_eq!(
labeled.atoms()[0].prop("_CIPNeighborOrder"),
Some("[4,3,2,1]")
);
}
#[test]
fn ciplabeler_assign_writes_sp2bond_labels_stereo_atoms_and_neighbor_order_like_rdkit() {
let molecule = sp2_bond_molecule(BondStereo::Cis, (4, 5));
let labeled = assign_cip_labels(&molecule, 0).unwrap();
assert_eq!(labeled.prop("_CIPComputed"), Some("1"));
assert_eq!(labeled.bonds()[0].prop("_CIPCode"), Some("Z"));
assert_eq!(labeled.bonds()[0].prop("_CIPNeighborOrder"), Some("[4,5]"));
assert_eq!(labeled.bonds()[0].stereo(), BondStereo::Cis);
assert_eq!(
labeled.bonds()[0].stereo_atoms(),
Some([crate::AtomId::new(4), crate::AtomId::new(5)])
);
}
#[test]
fn ciplabeler_assign_writes_atropisomer_bond_labels_like_rdkit() {
let molecule = atropisomer_bond_molecule(BondStereo::AtropCcw, false);
let labeled = assign_cip_labels(&molecule, 0).unwrap();
assert_eq!(labeled.prop("_CIPComputed"), Some("1"));
assert_eq!(labeled.bonds()[0].prop("_CIPCode"), Some("M"));
assert_eq!(labeled.bonds()[0].prop("_CIPNeighborOrder"), Some("[4,5]"));
}
#[test]
fn ciplabeler_assign_selected_masks_only_label_selected_atoms_and_bonds_like_rdkit() {
let mut molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
molecule.topology_block_mut().atoms[0].set_prop("_CIPCode", "old");
molecule.topology_block_mut().atoms[0].set_prop("_CIPNeighborOrder", "[0]");
let atom_mask = vec![false; molecule.num_atoms()];
let bond_mask = vec![false; molecule.num_bonds()];
let labeled = assign_cip_labels_for_indices(&molecule, &atom_mask, &bond_mask, 0).unwrap();
assert_eq!(labeled.prop("_CIPComputed"), Some("1"));
assert_eq!(labeled.atoms()[0].prop("_CIPCode"), Some("old"));
assert_eq!(labeled.atoms()[0].prop("_CIPNeighborOrder"), Some("[0]"));
let mut atom_mask = vec![false; molecule.num_atoms()];
atom_mask[0] = true;
let labeled = assign_cip_labels_for_indices(&molecule, &atom_mask, &bond_mask, 0).unwrap();
assert_eq!(labeled.atoms()[0].prop("_CIPCode"), Some("S"));
assert_eq!(
labeled.atoms()[0].prop("_CIPNeighborOrder"),
Some("[4,3,2,1]")
);
}
#[test]
fn ciplabeler_assign_max_recursive_iterations_does_not_disable_constitutional_fast_pass_like_rdkit()
{
let molecule = tetrahedral_molecule(ChiralTag::TetrahedralCcw);
let labeled = assign_cip_labels(&molecule, 1).unwrap();
assert_eq!(labeled.prop("_CIPComputed"), Some("1"));
assert_eq!(labeled.atoms()[0].prop("_CIPCode"), Some("S"));
assert_eq!(
labeled.atoms()[0].prop("_CIPNeighborOrder"),
Some("[4,3,2,1]")
);
}
}