zachcp/sssr.clj

## sssr.clj
(ns sssr.core
  (:gen-class))

(def pentagon2 {\A #{\B \E}, \B #{\A \C}
                \C #{\B \D}, \D #{\C \E}
                \E #{\D \A} })

(def hexagon2 {\A #{\B \F}, \B #{\A \C}
               \C #{\B \D}, \D #{\C \E}
               \E #{\D \F}, \F #{\E \A} })

(def four_six2 {\A #{\B \H},    \B #{\A \C}
                \C #{\B \D \H}, \D #{\C \E}
                \E #{\D \F},    \F #{\E \G}
                \G #{\H \F},    \H #{\A \C \G} })


(def three_five2 {\A #{\B \F},    \B #{\A \C \F}
                 \C #{\B \D},    \D #{\C \E}
                 \E #{\D \F},    \F #{\E \A \B} })


(def four_six_plus2 {\A #{\B \H \I}, \B #{\A \C}
                     \C #{\B \D \H}, \D #{\C \E}
                     \E #{\D \F},    \F #{\E \G}
                     \G #{\H \F},    \H #{\A \C \G}
                     \I #{\A \J},    \J #{\I} })


(defn dequeue!
  ; deque from http:\\stackoverflow.com\questions\8938330\clojure-swap-atom-dequeuing?rq=1
  [queue]
  (loop []
    (let [q     @queue
          value (peek q)
          nq    (pop q)]
      (if (compare-and-set! queue q nq)
        value
        (recur)))))


;try to put transformations in the form of a reducing function that works on a transient
; atoms are not seable so i was running to trouble with them transients are having trouble too

(defn prune   [Molecule]
  "take a Molecule graph and remore all atoms lacking 2 connections

  singles are the atoms with fewere than two connections
  atompairs are a vector of tuples corresponding to bonds of singles
  remove connection will remove the connections to an atom and the atom using atompairs

  the reduction will pass the Molecule through each pruning event. once there are no
  singels detected, the Molecule is returned"
  (let [singles    (filter #(< (count (second %)) 2) Molecule)
        atompairs  (partition 2 (flatten
                                  (for [[atm bond] singles]
                                    (for [b bond] [atm b]))))
        remove-connections (fn [mol keys12]
                             (let [[key1 key2] keys12
                                   nobond (update-in mol [key2] disj key1)
                                   noatom (dissoc nobond key1)]
                               noatom ))]
    (if (empty? singles)
      Molecule
      (let [pruned (reduce remove-connections Molecule atompairs)]
        (prune pruned)))))

(defn getring
  "breadth-first search for smallest ring modified BFS algorithm after Figueras
   all mutable state is kept in atoms:
      the atomqueue is used to kep track of which atom to evaluate
      the paths is a map of sets that keep track of all the paths starting at the rootnode
      the visited atom keeps track of atoms you've seen
      the ring is set to nil and awaits discovery

  Evaluation of smallest ring is assessed by merging your paths and making sure the intersection
  is only one."
  [Molecule rootnode]
  (let [ init-paths (into {} (for [[k v]  Molecule] [k #{}]))
        paths (atom (update-in init-paths [rootnode] conj rootnode))
		    atomqueue (atom (conj clojure.lang.PersistentQueue/EMPTY rootnode))
		    visited (atom #{rootnode})
		    ring (atom nil) ]
		 (while (nil? @ring)
       (let [atomcode (dequeue! atomqueue)
            connectedatoms (->> (get Molecule atomcode)
                                (remove @visited))]
         (doseq [atm connectedatoms]
           ;check for rings and return if a ring is found
           (let [currentpath (get @paths atomcode)
                nextpath  (get @paths atm)
                intersect (clojure.set/intersection currentpath nextpath)]
             (cond
               (and (> (count @paths) 1) (= 1 (count intersect)))
                 (reset! ring (clojure.set/union currentpath nextpath))
               (= 0 (count nextpath))
                 (do
                  (swap! paths update-in [atm] clojure.set/union currentpath)
                  (swap! paths update-in [atm] conj atm)
                  (swap! atomqueue conj atm)))))
         (swap! visited conj atomcode)))
    @ring))


(defn getrings
  "this funtion returns the smallest set of smallest rings of a Molecule graph
  it uses the prune funciton to eliminate unconnected atoms and uses the getrings
  funciton to find the smallest ring for each node

  unlike the Figueras algorithim this does not break the chains after finding a node
  therefore this will do more searches. A possible speed-up if needed soule be to implement
  chain breaking although the speedup would probably only benefit large molecules."
	[Molecule]
  (let [mol (prune Molecule)
        getrings1 (fn [rings x]
                    (let [ring (getring Molecule x)]
                      (clojure.set/union rings #{ring})))]
    (reduce getrings1 #{} (keys mol))))


(comment
  (prune four_six_plus2)
  (getring (prune hexagon2) \A)
  (getrings four_six2)
  (getrings four_six_plus2)
  )

## sssr.py
from collections import deque

pentagon = [{"key":"A", "connections":["B","E"]},
   {"key":"B", "connections": ["A", "C"]},
   {"key":"C", "connections": ["B", "D"]},
   {"key":"D", "connections": ["C", "E"]},
   {"key":"E", "connections": ["D", "A"]}]

hexagon = [{"key":"A","connections":["B", "E"]},
   {"key":"B","connections":["A", "C"]},
   {"key":"C","connections":["B", "D"]},
   {"key":"D","connections":["C", "E"]},
   {"key":"E","connections":["D", "A"]}]

four_six = [
{"key":"A","connections":["B", "H"]},
{"key":"B","connections":["A", "C"]},
{"key":"C","connections":["B", "D", "H"]},
{"key":"D","connections":["C", "E"]},
{"key":"E","connections":["D", "F"]},
{"key":"F","connections":["E", "G"]},
{"key":"G","connections":["F", "H"]},
{"key":"H","connections":["A", "C", "G"]}
]

three_five = [
{"key":"A","connections":["B", "F"]},
{"key":"B","connections":["A", "C", "F"]},
{"key":"C","connections":["B", "D"]},
{"key":"D","connections":["C", "E"]},
{"key":"E","connections":["D", "F"]},
{"key":"F","connections":["E", "A", "B"]}
]

four_six_plus = [
{"key":"A","connections":["B", "H", "I"]},
{"key":"B","connections":["A", "C"]},
{"key":"C","connections":["B", "D", "H"]},
{"key":"D","connections":["C", "E"]},
{"key":"E","connections":["D", "F"]},
{"key":"F","connections":["E", "G"]},
{"key":"G","connections":["F", "H"]},
{"key":"H","connections":["A", "C", "G"]},
{"key":"I","connections":["A","J"]},
{"key":"J","connections":["I"]}
]


def getring(Molecule, rootnode):
    # connections is just the table of connections in the molecule
    connections = { mol['key']: mol['connections'] for mol in Molecule}
    # ring members must be a member of the ring
    assert len(connections[rootnode])>1
    # paths is a dictionary of paths starting form the start node. It is initialized
    # to zero and will be used to updated as each member of the queue is analyzed
    paths = {mol['key']: set() for mol in Molecule}

    # atomqueue will keep track of atoms that have been visited and need to be processed
    atomqueue = deque()

    # visited is the atom keys that have been seen already
    visited = set()

    #when a ring is found, return the ring
    ring = None

    # initialize the root node
    atomqueue.append(rootnode)
    paths[rootnode].add(rootnode)
    visited.add(rootnode)

    while not ring:
        print "atomqueue: {}".format(atomqueue)
        print "visited: {}".format(visited)
        print "paths: {}".format(paths)

        try:
            atomcode = atomqueue.popleft()
            #print "atomcode: {}".format(atomcode)

            connectedatoms = [a for a in connections[atomcode] if a not in visited]
            #print "atoms connected: {}".format(connectedatoms)
            for atom in connectedatoms:
                print "atom: {}".format(atom)
                #check for rings and return if a ring is found
                currentpath = paths[atomcode]
                nextpath  = paths[atom]
                intersect = currentpath.intersection(nextpath)
                #print "current path: {}".format(currentpath)
                #print "next path: {}".format(currentpath)
                if len(nextpath) > 0 and len(intersect) == 1:
                    return currentpath.union(nextpath)

                # update path only if the target path is empty
                # to avoid incorporating longer paths
                if len(nextpath) == 0:
                    paths[atom] = currentpath.copy()
                    paths[atom].add(atom)

                #add atoms to queue
                atomqueue.append(atom)

            #update visited
            visited.add(atomcode)
        except:
            raise Error("Problem Here. You should never get here!")

def getrings(Molecule):
    """
    find atoms that are connected by only a single bond,
    eliminate them and the connections to them by other atoms until
    there are only atoms that have at least two connections
    then find the rings on them
    """

    # copy the Molecule for local editing
    tempmol = list(Molecule)
    print tempmol
    # find the atoms that are initially less than 2 bond connections
    singles_index = [i for i,a in enumerate(tempmol) if len(a['connections']) < 2]


    # eliminate all atoms that are not part of a ring
    while len(singles_index) > 0:
        for idx in singles_index:
            key = tempmol[idx]['key']
            connections = tempmol[idx]['connections']

            # remove connections
            for con in connections:
                #print "index: {}\nkey:{}\nconnection:{}".format(idx,key,con)
                index = [i for i,x in enumerate(tempmol) if x['key'] == con][0]
                tempmol[index]['connections'].remove(key)

            # remove atom
            tempmol.pop(idx)
            print tempmol

            # reset singles index
            singles_index = [i for i,a in enumerate(tempmol) if len(a['connections']) < 2]
            print len(singles_index)

    # calculate ring on the trimmed molceule
    rings = set()
    for atom in tempmol:
        if len(atom['connections']) > 1:
            ring = frozenset(getring(Molecule, atom["key"]))
            print atom, ring
            if ring not in rings:
                rings.add(ring)
    return rings


#getring(pentagon, 'B') #expect six
#getring(hexagon, 'B')  # expect five
#getring(four_six, 'A') #expect four yes
#getring(four_six, 'B') #expect four yes
#getring(four_six, 'C') #expect four yes
#getring(four_six, 'D') #expect six yes
#getring(four_six, 'E') #expect six yes
#getring(four_six, 'F') #expect six yes
#getring(four_six, 'G') #expect siz yes
#getring(four_six, 'H') #expect four yes
#getring(four_six_plus, 'I') #expect assertion for length
#getrings(four_six) # correctly gets rings
#getrings(three_five) # incorrect for C,E,F
#getrings(four_six_plus) # should trim atoms and correctly get rings
	(ns sssr.core
	(:gen-class))

	(def pentagon2 {\A #{\B \E}, \B #{\A \C}
	\C #{\B \D}, \D #{\C \E}
	\E #{\D \A} })

	(def hexagon2 {\A #{\B \F}, \B #{\A \C}
	\C #{\B \D}, \D #{\C \E}
	\E #{\D \F}, \F #{\E \A} })

	(def four_six2 {\A #{\B \H}, \B #{\A \C}
	\C #{\B \D \H}, \D #{\C \E}
	\E #{\D \F}, \F #{\E \G}
	\G #{\H \F}, \H #{\A \C \G} })


	(def three_five2 {\A #{\B \F}, \B #{\A \C \F}
	\C #{\B \D}, \D #{\C \E}
	\E #{\D \F}, \F #{\E \A \B} })


	(def four_six_plus2 {\A #{\B \H \I}, \B #{\A \C}
	\C #{\B \D \H}, \D #{\C \E}
	\E #{\D \F}, \F #{\E \G}
	\G #{\H \F}, \H #{\A \C \G}
	\I #{\A \J}, \J #{\I} })


	(defn dequeue!
	; deque from http:\\stackoverflow.com\questions\8938330\clojure-swap-atom-dequeuing?rq=1
	[queue]
	(loop []
	(let [q @queue
	value (peek q)
	nq (pop q)]
	(if (compare-and-set! queue q nq)
	value
	(recur)))))


	;try to put transformations in the form of a reducing function that works on a transient
	; atoms are not seable so i was running to trouble with them transients are having trouble too

	(defn prune [Molecule]
	"take a Molecule graph and remore all atoms lacking 2 connections

	singles are the atoms with fewere than two connections
	atompairs are a vector of tuples corresponding to bonds of singles
	remove connection will remove the connections to an atom and the atom using atompairs

	the reduction will pass the Molecule through each pruning event. once there are no
	singels detected, the Molecule is returned"
	(let [singles (filter #(< (count (second %)) 2) Molecule)
	atompairs (partition 2 (flatten
	(for [[atm bond] singles]
	(for [b bond] [atm b]))))
	remove-connections (fn [mol keys12]
	(let [[key1 key2] keys12
	nobond (update-in mol [key2] disj key1)
	noatom (dissoc nobond key1)]
	noatom ))]
	(if (empty? singles)
	Molecule
	(let [pruned (reduce remove-connections Molecule atompairs)]
	(prune pruned)))))

	(defn getring
	"breadth-first search for smallest ring modified BFS algorithm after Figueras
	all mutable state is kept in atoms:
	the atomqueue is used to kep track of which atom to evaluate
	the paths is a map of sets that keep track of all the paths starting at the rootnode
	the visited atom keeps track of atoms you've seen
	the ring is set to nil and awaits discovery

	Evaluation of smallest ring is assessed by merging your paths and making sure the intersection
	is only one."
	[Molecule rootnode]
	(let [ init-paths (into {} (for [[k v] Molecule] [k #{}]))
	paths (atom (update-in init-paths [rootnode] conj rootnode))
	atomqueue (atom (conj clojure.lang.PersistentQueue/EMPTY rootnode))
	visited (atom #{rootnode})
	ring (atom nil) ]
	(while (nil? @ring)
	(let [atomcode (dequeue! atomqueue)
	connectedatoms (->> (get Molecule atomcode)
	(remove @visited))]
	(doseq [atm connectedatoms]
	;check for rings and return if a ring is found
	(let [currentpath (get @paths atomcode)
	nextpath (get @paths atm)
	intersect (clojure.set/intersection currentpath nextpath)]
	(cond
	(and (> (count @paths) 1) (= 1 (count intersect)))
	(reset! ring (clojure.set/union currentpath nextpath))
	(= 0 (count nextpath))
	(do
	(swap! paths update-in [atm] clojure.set/union currentpath)
	(swap! paths update-in [atm] conj atm)
	(swap! atomqueue conj atm)))))
	(swap! visited conj atomcode)))
	@ring))


	(defn getrings
	"this funtion returns the smallest set of smallest rings of a Molecule graph
	it uses the prune funciton to eliminate unconnected atoms and uses the getrings
	funciton to find the smallest ring for each node

	unlike the Figueras algorithim this does not break the chains after finding a node
	therefore this will do more searches. A possible speed-up if needed soule be to implement
	chain breaking although the speedup would probably only benefit large molecules."
	[Molecule]
	(let [mol (prune Molecule)
	getrings1 (fn [rings x]
	(let [ring (getring Molecule x)]
	(clojure.set/union rings #{ring})))]
	(reduce getrings1 #{} (keys mol))))


	(comment
	(prune four_six_plus2)
	(getring (prune hexagon2) \A)
	(getrings four_six2)
	(getrings four_six_plus2)
	)
	from collections import deque

	pentagon = [{"key":"A", "connections":["B","E"]},
	{"key":"B", "connections": ["A", "C"]},
	{"key":"C", "connections": ["B", "D"]},
	{"key":"D", "connections": ["C", "E"]},
	{"key":"E", "connections": ["D", "A"]}]

	hexagon = [{"key":"A","connections":["B", "E"]},
	{"key":"B","connections":["A", "C"]},
	{"key":"C","connections":["B", "D"]},
	{"key":"D","connections":["C", "E"]},
	{"key":"E","connections":["D", "A"]}]

	four_six = [
	{"key":"A","connections":["B", "H"]},
	{"key":"B","connections":["A", "C"]},
	{"key":"C","connections":["B", "D", "H"]},
	{"key":"D","connections":["C", "E"]},
	{"key":"E","connections":["D", "F"]},
	{"key":"F","connections":["E", "G"]},
	{"key":"G","connections":["F", "H"]},
	{"key":"H","connections":["A", "C", "G"]}
	]

	three_five = [
	{"key":"A","connections":["B", "F"]},
	{"key":"B","connections":["A", "C", "F"]},
	{"key":"C","connections":["B", "D"]},
	{"key":"D","connections":["C", "E"]},
	{"key":"E","connections":["D", "F"]},
	{"key":"F","connections":["E", "A", "B"]}
	]

	four_six_plus = [
	{"key":"A","connections":["B", "H", "I"]},
	{"key":"B","connections":["A", "C"]},
	{"key":"C","connections":["B", "D", "H"]},
	{"key":"D","connections":["C", "E"]},
	{"key":"E","connections":["D", "F"]},
	{"key":"F","connections":["E", "G"]},
	{"key":"G","connections":["F", "H"]},
	{"key":"H","connections":["A", "C", "G"]},
	{"key":"I","connections":["A","J"]},
	{"key":"J","connections":["I"]}
	]


	def getring(Molecule, rootnode):
	# connections is just the table of connections in the molecule
	connections = { mol['key']: mol['connections'] for mol in Molecule}
	# ring members must be a member of the ring
	assert len(connections[rootnode])>1
	# paths is a dictionary of paths starting form the start node. It is initialized
	# to zero and will be used to updated as each member of the queue is analyzed
	paths = {mol['key']: set() for mol in Molecule}

	# atomqueue will keep track of atoms that have been visited and need to be processed
	atomqueue = deque()

	# visited is the atom keys that have been seen already
	visited = set()

	#when a ring is found, return the ring
	ring = None

	# initialize the root node
	atomqueue.append(rootnode)
	paths[rootnode].add(rootnode)
	visited.add(rootnode)

	while not ring:
	print "atomqueue: {}".format(atomqueue)
	print "visited: {}".format(visited)
	print "paths: {}".format(paths)

	try:
	atomcode = atomqueue.popleft()
	#print "atomcode: {}".format(atomcode)

	connectedatoms = [a for a in connections[atomcode] if a not in visited]
	#print "atoms connected: {}".format(connectedatoms)
	for atom in connectedatoms:
	print "atom: {}".format(atom)
	#check for rings and return if a ring is found
	currentpath = paths[atomcode]
	nextpath = paths[atom]
	intersect = currentpath.intersection(nextpath)
	#print "current path: {}".format(currentpath)
	#print "next path: {}".format(currentpath)
	if len(nextpath) > 0 and len(intersect) == 1:
	return currentpath.union(nextpath)

	# update path only if the target path is empty
	# to avoid incorporating longer paths
	if len(nextpath) == 0:
	paths[atom] = currentpath.copy()
	paths[atom].add(atom)

	#add atoms to queue
	atomqueue.append(atom)

	#update visited
	visited.add(atomcode)
	except:
	raise Error("Problem Here. You should never get here!")

	def getrings(Molecule):
	"""
	find atoms that are connected by only a single bond,
	eliminate them and the connections to them by other atoms until
	there are only atoms that have at least two connections
	then find the rings on them
	"""

	# copy the Molecule for local editing
	tempmol = list(Molecule)
	print tempmol
	# find the atoms that are initially less than 2 bond connections
	singles_index = [i for i,a in enumerate(tempmol) if len(a['connections']) < 2]


	# eliminate all atoms that are not part of a ring
	while len(singles_index) > 0:
	for idx in singles_index:
	key = tempmol[idx]['key']
	connections = tempmol[idx]['connections']

	# remove connections
	for con in connections:
	#print "index: {}\nkey:{}\nconnection:{}".format(idx,key,con)
	index = [i for i,x in enumerate(tempmol) if x['key'] == con][0]
	tempmol[index]['connections'].remove(key)

	# remove atom
	tempmol.pop(idx)
	print tempmol

	# reset singles index
	singles_index = [i for i,a in enumerate(tempmol) if len(a['connections']) < 2]
	print len(singles_index)

	# calculate ring on the trimmed molceule
	rings = set()
	for atom in tempmol:
	if len(atom['connections']) > 1:
	ring = frozenset(getring(Molecule, atom["key"]))
	print atom, ring
	if ring not in rings:
	rings.add(ring)
	return rings


	#getring(pentagon, 'B') #expect six
	#getring(hexagon, 'B') # expect five
	#getring(four_six, 'A') #expect four yes
	#getring(four_six, 'B') #expect four yes
	#getring(four_six, 'C') #expect four yes
	#getring(four_six, 'D') #expect six yes
	#getring(four_six, 'E') #expect six yes
	#getring(four_six, 'F') #expect six yes
	#getring(four_six, 'G') #expect siz yes
	#getring(four_six, 'H') #expect four yes
	#getring(four_six_plus, 'I') #expect assertion for length
	#getrings(four_six) # correctly gets rings
	#getrings(three_five) # incorrect for C,E,F
	#getrings(four_six_plus) # should trim atoms and correctly get rings