1 package Graph::CyclesDetection;
   2 #
   3 # File: CyclesDetection.pm
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2024 Manish Sud. All rights reserved.
   7 #
   8 # This file is part of MayaChemTools.
   9 #
  10 # MayaChemTools is free software; you can redistribute it and/or modify it under
  11 # the terms of the GNU Lesser General Public License as published by the Free
  12 # Software Foundation; either version 3 of the License, or (at your option) any
  13 # later version.
  14 #
  15 # MayaChemTools is distributed in the hope that it will be useful, but without
  16 # any warranty; without even the implied warranty of merchantability of fitness
  17 # for a particular purpose.  See the GNU Lesser General Public License for more
  18 # details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public License
  21 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  22 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  23 # Boston, MA, 02111-1307, USA.
  24 #
  25 
  26 use strict;
  27 use Carp;
  28 use Exporter;
  29 use Graph;
  30 use Graph::Path;
  31 use Graph::PathGraph;
  32 use BitVector;
  33 
  34 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  35 
  36 @ISA = qw(Exporter);
  37 @EXPORT = qw();
  38 @EXPORT_OK = qw();
  39 
  40 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
  41 
  42 # Setup class variables...
  43 my($ClassName);
  44 _InitializeClass();
  45 
  46 # Overload Perl functions...
  47 use overload '""' => 'StringifyCyclesDetection';
  48 
  49 # Class constructor...
  50 sub new {
  51   my($Class, $Graph) = @_;
  52 
  53   # Initialize object...
  54   my $This = {};
  55   bless $This, ref($Class) || $Class;
  56   $This->_InitializeCyclesDetection($Graph);
  57 
  58   return $This;
  59 }
  60 
  61 # Initialize object data...
  62 sub _InitializeCyclesDetection {
  63   my($This, $Graph) = @_;
  64 
  65   $This->{Graph} = $Graph;
  66 
  67   # Cycles list...
  68   @{$This->{AllCyclicPaths}} = ();
  69 
  70   # Cyclic paths which are not part of any other cycle...
  71   @{$This->{IndependentCyclicPaths}} = ();
  72 
  73   return $This;
  74 }
  75 
  76 # Initialize class ...
  77 sub _InitializeClass {
  78   #Class name...
  79   $ClassName = __PACKAGE__;
  80 }
  81 
  82 # Detect all cycles in graph...
  83 #
  84 sub DetectCycles {
  85   my($This) = @_;
  86 
  87   return $This->DetectCyclesUsingCollapsingPathGraphMethodology();
  88 }
  89 
  90 # Detect all cycles in the graph using collapsing path graph [Ref 31] methodology...
  91 #
  92 # Note:
  93 #   . For topologically complex graphs containing large number of cycles,
  94 #     CollapseVertexAndCollectCyclicPathsDetectCycles method implemented in
  95 #     PathGraph can time out in which case no cycles are detected or
  96 #     assigned.
  97 #
  98 sub DetectCyclesUsingCollapsingPathGraphMethodology {
  99   my($This) = @_;
 100   my($PathGraph);
 101 
 102   # Create a path graph...
 103   $PathGraph = new Graph::PathGraph($This->{Graph});
 104 
 105   # For a vertex to be in a cycle, its degree must be >=2. So delete vertices recursively
 106   # till all vertices with degree less than 2 have been deleted...
 107   $PathGraph->DeleteVerticesWithDegreeLessThan(2);
 108 
 109   # Setup a VertexID and EdgeID to index map usage during retrieval of independent cycles to
 110   # avoid going over all vertices in all cylic paths later...
 111   #
 112   my($VertexIDsToIndicesRef, $LargestVertexIndex, $EdgeIDsToIndicesRef, $LargestEdgeIDIndex);
 113   ($VertexIDsToIndicesRef, $LargestVertexIndex) = $This->_SetupVertexIDsToIndicesMap($PathGraph);
 114   ($EdgeIDsToIndicesRef, $LargestEdgeIDIndex) = $This->_SetupEdgeIDsToIndicesMap($PathGraph);
 115 
 116   # Recursively collapse vertices with lowest degree...
 117   my($VertexID, $CycleVertexID);
 118   while ($VertexID = $PathGraph->GetVertexWithSmallestDegree()) {
 119       if (!$PathGraph->CollapseVertexAndCollectCyclicPaths($VertexID)) {
 120         # Cycles detection didn't finish...
 121         return undef;
 122       }
 123   }
 124 
 125   # Get detected cycles and save 'em sorted by size...
 126   push @{$This->{AllCyclicPaths}}, sort { $a->GetLength() <=> $b->GetLength() } $PathGraph->GetCyclicPaths();
 127 
 128   # Retrieve independent cyclic paths...
 129   return $This->_RetrieveIndependentCycles($VertexIDsToIndicesRef, $LargestVertexIndex, $EdgeIDsToIndicesRef, $LargestEdgeIDIndex);
 130 }
 131 
 132 # Retrieve and save independent cyclic paths...
 133 #
 134 # Set of independent cycles identified using this method doesn't correspond to basis set of
 135 # rings or smallest set of smallest rings (SSSR) [ Refs 29-30 ]; instead, set of cycles identified
 136 # as independent cycles simply correspond to cycles which contain no other cycle as their
 137 # subcycles and can't be described as linear combination of smaller cycles. And it also happen
 138 # to contain all the rings in basis set of rings and SSSR. In other words, it's a superset of basis set
 139 # of cycles and SSSR. For example, six four membered cycles are identified for cubane which is one
 140 # more than the basis set of cycles.
 141 #
 142 sub _RetrieveIndependentCycles {
 143   my($This, $VertexIDsToIndicesRef, $LargestVertexIndex, $EdgeIDsToIndicesRef, $LargestEdgeIDIndex) = @_;
 144 
 145   # Only 1 or 0 cyclic paths...
 146   if (@{$This->{AllCyclicPaths}} <= 1) {
 147     push @{$This->{IndependentCyclicPaths}}, @{$This->{AllCyclicPaths}};
 148     return $This;
 149   }
 150 
 151   # Setup bit vectors for each cyclic path with size of each bit vector corresponding to
 152   # maximum vertex index plus one...
 153   my($CyclicPath, $BitVector, @BitNums, @CyclicPathBitVectors, @CyclicPathEdgeIDsBitVectors);
 154 
 155   @CyclicPathBitVectors = (); @CyclicPathEdgeIDsBitVectors = ();
 156 
 157   # Set bits corresponding to VertexIDs EdgeIDs using their indices...
 158   for $CyclicPath (@{$This->{AllCyclicPaths}}) {
 159     $BitVector = new BitVector($LargestVertexIndex);
 160     @BitNums = map { $VertexIDsToIndicesRef->{$_} } $CyclicPath->GetVertices();
 161     $BitVector->SetBits(@BitNums);
 162     push @CyclicPathBitVectors, $BitVector;
 163 
 164     $BitVector = new BitVector($LargestEdgeIDIndex);
 165     @BitNums = map { $EdgeIDsToIndicesRef->{$_} } $This->_GetPathEdgeIDs($CyclicPath);
 166     $BitVector->SetBits(@BitNums);
 167     push @CyclicPathEdgeIDsBitVectors, $BitVector;
 168   }
 169 
 170   # First smallest cyclic path always ends up as an independent cyclic path...
 171   push @{$This->{IndependentCyclicPaths}}, $This->{AllCyclicPaths}[0];
 172 
 173   # Retrieve other independent cyclic paths...
 174   my($CurrentIndex, $PreviousIndex, $CurrentCyclicPath, $PreviousCyclicPath, $CurrentPathLength, $PreviousPathLength, $CurrentBitVector, $PreviousBitVector, $CurrentAndPreviousBitVectpor, $AllPreviousSmallerPathsBitVector, $AllPreviousSmallerPathsEdgeIDsBitVector, $CurrentEdgeIDsBitVector, $AndBitVector, %SmallerPathAlreadyAdded, %SkipPath);
 175 
 176   %SkipPath = ();
 177   %SmallerPathAlreadyAdded = ();
 178   $AllPreviousSmallerPathsBitVector = new BitVector($LargestVertexIndex);
 179   $AllPreviousSmallerPathsEdgeIDsBitVector = new BitVector($LargestEdgeIDIndex);
 180 
 181   CURRENT: for $CurrentIndex (1 .. $#{$This->{AllCyclicPaths}}) {
 182     if (exists $SkipPath{$CurrentIndex}) {
 183       next CURRENT;
 184     }
 185     $CurrentCyclicPath = $This->{AllCyclicPaths}[$CurrentIndex];
 186     $CurrentBitVector = $CyclicPathBitVectors[$CurrentIndex];
 187     $CurrentPathLength = $CurrentCyclicPath->GetLength();
 188 
 189     PREVIOUS: for $PreviousIndex (0 .. ($CurrentIndex - 1)) {
 190       if (exists $SkipPath{$PreviousIndex}) {
 191         next PREVIOUS;
 192       }
 193       $PreviousCyclicPath = $This->{AllCyclicPaths}[$PreviousIndex];
 194       $PreviousBitVector = $CyclicPathBitVectors[$PreviousIndex];
 195 
 196       # Is previous path a subset of current path?
 197       $CurrentAndPreviousBitVectpor = $PreviousBitVector &  $CurrentBitVector;
 198       if ($PreviousBitVector->GetNumOfSetBits() == $CurrentAndPreviousBitVectpor->GetNumOfSetBits()) {
 199         # Current path doesn't qualify an independent path...
 200         $SkipPath{$CurrentIndex} = 1;
 201         next CURRENT;
 202       }
 203 
 204       $PreviousPathLength = $PreviousCyclicPath->GetLength();
 205       if ($PreviousPathLength < $CurrentPathLength) {
 206         # Mark cycle vertices seen in cyclic paths with length smaller than current path...
 207         if (! exists $SmallerPathAlreadyAdded{$PreviousIndex}) {
 208           $SmallerPathAlreadyAdded{$PreviousIndex} = 1;
 209           $AllPreviousSmallerPathsBitVector = $AllPreviousSmallerPathsBitVector | $PreviousBitVector;
 210           $AllPreviousSmallerPathsEdgeIDsBitVector = $AllPreviousSmallerPathsEdgeIDsBitVector | $CyclicPathEdgeIDsBitVectors[$PreviousIndex];
 211         }
 212       }
 213     }
 214     if ($AllPreviousSmallerPathsBitVector->GetNumOfSetBits()) {
 215       # Is current path a linear combination of smaller paths?
 216       $AndBitVector = $AllPreviousSmallerPathsBitVector &  $CurrentBitVector;
 217       if ($CurrentBitVector->GetNumOfSetBits() == $AndBitVector->GetNumOfSetBits()) {
 218         # Are all edges in the current path already present in smaller paths?
 219         $CurrentEdgeIDsBitVector = $CyclicPathEdgeIDsBitVectors[$CurrentIndex];
 220         $AndBitVector = $AllPreviousSmallerPathsEdgeIDsBitVector &  $CurrentEdgeIDsBitVector;
 221 
 222         if ($CurrentEdgeIDsBitVector->GetNumOfSetBits() == $AndBitVector->GetNumOfSetBits()) {
 223           # Current path doesn't qualify an independent path...
 224           $SkipPath{$CurrentIndex} = 1;
 225           next CURRENT;
 226         }
 227       }
 228     }
 229     # It's an independent cyclic path...
 230     push @{$This->{IndependentCyclicPaths}}, $CurrentCyclicPath;
 231   }
 232   return $This;
 233 }
 234 
 235 # Setup a VertexID to index map...
 236 #
 237 sub _SetupVertexIDsToIndicesMap {
 238   my($This, $PathGraph) = @_;
 239   my($LargestVertexIndex, @VertexIDs, %VertexIDsMap);
 240 
 241   %VertexIDsMap = (); @VertexIDs = (); $LargestVertexIndex = 0;
 242 
 243   @VertexIDs = $PathGraph->GetVertices();
 244   if (! @VertexIDs) {
 245     return (\%VertexIDsMap, $LargestVertexIndex);
 246   }
 247   @VertexIDsMap{ @VertexIDs } = (0 .. $#VertexIDs);
 248   $LargestVertexIndex = scalar @VertexIDs;
 249 
 250   return (\%VertexIDsMap, $LargestVertexIndex);
 251 }
 252 
 253 # Setup a Edge to index map using paths associated to egdes in an intial
 254 # path graph...
 255 #
 256 sub _SetupEdgeIDsToIndicesMap {
 257   my($This, $PathGraph) = @_;
 258   my($Path, $LargestEdgeIndex, @EdgeIDs, %EdgeIDsMap);
 259 
 260   %EdgeIDsMap = (); @EdgeIDs = (); $LargestEdgeIndex = 0;
 261 
 262   @EdgeIDs = ();
 263   for $Path ($PathGraph->GetPaths()) {
 264     push @EdgeIDs, $This->_GetPathEdgeIDs($Path);
 265   }
 266 
 267   if (! @EdgeIDs) {
 268     return (\%EdgeIDsMap, $LargestEdgeIndex);
 269   }
 270 
 271   @EdgeIDsMap{ @EdgeIDs } = (0 .. $#EdgeIDs);
 272   $LargestEdgeIndex = scalar @EdgeIDs;
 273 
 274   return (\%EdgeIDsMap, $LargestEdgeIndex);
 275 }
 276 
 277 # Get path edge IDs or number of edges. The edge IDs are generated from
 278 # edge vertices and correpond to VertexID1-VertexID2 where VertexID1 <=
 279 # VertexID2.
 280 #
 281 sub _GetPathEdgeIDs {
 282   my($This, $Path) = @_;
 283   my(@EdgesVertexIDs, @EdgeIDs);
 284 
 285   @EdgesVertexIDs = (); @EdgeIDs = ();
 286   @EdgesVertexIDs = $Path->GetEdges();
 287   if (!@EdgesVertexIDs) {
 288     return wantarray ? @EdgeIDs : (scalar @EdgeIDs);
 289   }
 290 
 291   # Set up edge IDs...
 292   my($Index, $VertexID1, $VertexID2, $EdgeID);
 293 
 294   for ($Index = 0; $Index < $#EdgesVertexIDs; $Index += 2) {
 295     $VertexID1 = $EdgesVertexIDs[$Index]; $VertexID2 = $EdgesVertexIDs[$Index + 1];
 296     $EdgeID = ($VertexID1 <= $VertexID2) ? "$VertexID1-$VertexID2" : "$VertexID2-$VertexID1";
 297     push @EdgeIDs, $EdgeID;
 298   }
 299 
 300   return wantarray ? @EdgeIDs : (scalar @EdgeIDs);
 301 }
 302 
 303 # Return an array containing references to cyclic paths or number of cylic paths...
 304 #
 305 sub GetAllCyclicPaths {
 306   my($This) = @_;
 307 
 308   return wantarray ? @{$This->{AllCyclicPaths}} : scalar @{$This->{AllCyclicPaths}};
 309 }
 310 
 311 # Get cyclic paths which are independent. In otherwords, cycles which don't contain any other
 312 # cycle as their subset...
 313 #
 314 sub GetIndependentCyclicPaths {
 315   my($This) = @_;
 316 
 317   return wantarray ? @{$This->{IndependentCyclicPaths}} : scalar @{$This->{IndependentCyclicPaths}};
 318 }
 319 
 320 # Return a string containg data for CyclesDetection object...
 321 sub StringifyCyclesDetection {
 322   my($This) = @_;
 323   my($CyclesDetectionString, $CyclesCount, $CyclicPath);
 324 
 325   $CyclesCount = @{$This->{AllCyclicPaths}};
 326   $CyclesDetectionString = "AllCycles: Count - $CyclesCount";
 327   for $CyclicPath (@{$This->{AllCyclicPaths}}) {
 328     $CyclesDetectionString .= "; Cycle: " . join('-', $CyclicPath->GetVertices());
 329   }
 330 
 331   $CyclesCount = @{$This->{IndependentCyclicPaths}};
 332   $CyclesDetectionString .= "\nIndependentCycles: Count - $CyclesCount";
 333   for $CyclicPath (@{$This->{IndependentCyclicPaths}}) {
 334     $CyclesDetectionString .= "; Cycle: " . join('-', $CyclicPath->GetVertices());
 335   }
 336 
 337   return $CyclesDetectionString;
 338 }
 339 
 340 # Return a reference to new cycle detection object...
 341 sub Copy {
 342   my($This) = @_;
 343   my($NewCyclesDetection);
 344 
 345   $NewCyclesDetection = Storable::dclone($This);
 346 
 347   return $NewCyclesDetection;
 348 }
 349