Fix minimizer and sketch issues

src/chem/minimizer.cc

@@ -542,8 +542,24 @@ double 	Minimizer_O::dbrent(	double ax, double bx, double cx,
     if (fabs(_e)>tol1){
       d1=2.0*(_b-_a);
       d2=d1;
-      if (dw != dx ) d1=(w-x)*dx/(dx-dw); // Secant method, first on one, then on
-      if (dv != dx ) d2=(v-x)*dx/(dx-dv); // the other point
+
+      double numerator;
+      double denomenator;
+      double limit;
+      //if (dw != dx ) d1=(w-x)*dx/(dx-dw); // Original code Secant method, first on one, then on
+#if 0
+# define WARN_BAD_VAL(xx) if (std::isnan(xx) || std::isinf(xx)) printf("%s:%d:%s problem with " #xx " - it is %e\n", __FILE__, __LINE__, __FUNCTION__, xx );
+#else
+# define WARN_BAD_VAL(xx)
+#endif
+      WARN_BAD_VAL(dw);
+      WARN_BAD_VAL(dx);
+      WARN_BAD_VAL(dv);
+      if ((dw-dx) != 0.0) d1=(w-x)*dx/(dx-dw); // Secant method, first on one, then on
+
+      // if (dv!=dx) d2=(v-x)*dx/(dx-dv); // the other point
+      if ((dv-dx) != 0.0) d2=(v-x)*dx/(dx-dv); // the other point
+
 	    // Which of these two estimates of d shall we take? We will insist that
 	    // they are within the bracket, and on the side pointed to by the
 	    // derivative at x

src/lisp/cando/conditions.lisp

@@ -26,9 +26,9 @@
 (in-package :chem)
 
 (defun linear-atoms-reporter (condition stream)
-  (let ((indices (chem::indices condition))
+  (let ((indexes (chem::indexes condition))
         (coordinates (coordinates condition)))
-    (let ((positions (loop for index in indices
+    (let ((positions (loop for index in indexes
                            for xpos = (elt coordinates index)
                            for ypos = (elt coordinates (+ 1 index))
                            for zpos = (elt coordinates (+ 2 index))
@@ -48,7 +48,7 @@
   ((atoms :initarg :atoms :reader chem:atoms)
    (coordinates :initarg :coordinates :reader coordinates)
    (simd-index :initarg :simd-index :reader simd-index)
-   (indices :initarg :indices :reader chem::indices))
+   (indexes :initarg :indexes :reader chem::indexes))
   (:report linear-atoms-reporter))
 
 (define-condition chem:linear-angle-error (linear-atoms-error) ())

src/lisp/cando/molecules.lisp

@@ -191,7 +191,8 @@ Example:  (set-stereoisomer-mapping *agg* '((:C1 :R) (:C2 :S))"
   (restart-case
       (with-handle-linear-angles-dihedrals (:verbose verbose)
         (multiple-value-bind (pos total-energy)
-            (chem:minimize minimizer)
+            (ext:with-float-traps-masked (:underflow :overflow :invalid :inexact :divide-by-zero)
+              (chem:minimize minimizer))
           (unless total-energy
             (break "total-energy was NIL"))
           (values pos total-energy)

src/lisp/sketch2d/render-svg.lisp

@@ -1,5 +1,9 @@
 (in-package :sketch2d)
 
+
+
+(defparameter *angles* '(0 180 90 270 45 135 215 305))
+(defparameter *multiplier* 1.5)
 (defparameter *number-of-hashes-in-hashed-bond* 5)
 (defparameter *show-all* nil)
 (defparameter *show-names* nil)
@@ -61,6 +65,7 @@
 
 (defclass sketch-svg ()
   ((molecule :initarg :molecule :accessor molecule)
+   (optimized-annotations :initform nil :initarg :optimized-annotations :reader optimized-annotations)
    (sketch2d* :initarg :sketch2d* :accessor sketch2d*)
    (height :initarg :height :accessor height)
    (width :initarg :width :accessor width)
@@ -191,6 +196,57 @@ This will place the calculated bond on one or the other side of the x1,y1-x2,y2
           (setf (gethash atom atom-ht) new-atom-node)
           new-atom-node))))
 
+(defclass optimized-annotation (annotation)
+  (
+   (average-pos :initarg :average-pos :reader average-pos)
+   (closest-annotation-pos :initarg :closest-annotation-pos :reader closest-annotation-pos)
+   ))
+
+(defun optimize-annotations (sketch2d atom-nodes atoms-to-nodes)
+  (let* ((annotations (annotations sketch2d))
+         (molecule (molecule sketch2d))
+         (center-pos (let ((sum-pos (geom:vec 0.0 0.0 0.0)))
+                       (loop for atom-node in atom-nodes
+                             do (setf sum-pos (geom:v+ sum-pos (pos atom-node))))
+                       (geom:v* sum-pos (/ 1.0 (length atom-nodes)))))
+         (radius (let ((max-dist most-negative-single-float))
+                   (loop for atom-node in atom-nodes
+                         for dist = (geom:vlength (geom:v- (pos atom-node) center-pos))
+                         do (setf max-dist (max dist max-dist)))
+                   max-dist)))
+    (loop with used-angles = (make-array (length *angles*) :element-type 'bit :initial-element 0)
+          for annotation in annotations
+          for average-pos = (let ((pos-sum (geom:vec 0.0 0.0 0.0)))
+                              (loop for atom-name in (atom-names annotation)
+                                    for atm = (chem:first-atom-with-name molecule atom-name)
+                                    for atom-node = (gethash atm atoms-to-nodes)
+                                    for pos = (pos atom-node)
+                                    do (setf pos-sum (geom:v+ pos pos-sum)))
+                              (let ((pp (geom:v* pos-sum (/ 1.0 (length (atom-names annotation))))))
+                                (cons (geom:vx pp) (geom:vy pp))))
+          for closest-annotation-pos = (let ((min-dist most-positive-single-float)
+                                             best-x best-y best-angle-index)
+                                         (loop for angle-index from 0
+                                               for angle-deg in *angles*
+                                               for angle-rad = (* 0.0174533 angle-deg)
+                                               for y-pos = (+ (geom:vy center-pos) (* *multiplier* radius (sin angle-rad)))
+                                               for x-pos = (+ (geom:vx center-pos) (* *multiplier* radius (cos angle-rad)))
+                                               for dist = (sqrt (+
+                                                                 (expt (- y-pos (cdr average-pos)) 2)
+                                                                 (expt (- x-pos (car average-pos)) 2)))
+                                               when (and (= (aref used-angles angle-index) 0) (< dist min-dist))
+                                                 do (setf min-dist dist
+                                                          best-angle-index angle-index
+                                                          best-x x-pos
+                                                          best-y y-pos))
+                                         (setf (aref used-angles best-angle-index) 1)
+                                         (cons best-x best-y))
+          collect (make-instance 'optimized-annotation
+                                 :text (text annotation)
+                                 :atom-names (atom-names annotation)
+                                 :average-pos average-pos
+                                 :closest-annotation-pos closest-annotation-pos))))
+
 (defun generate-sketch (sketch2d width height)
   "Return the bonds and atoms that we need to render"
   (let ((molecule (molecule sketch2d))
@@ -442,7 +498,7 @@ This will place the calculated bond on one or the other side of the x1,y1-x2,y2
     (loop for bond-node in (bond-nodes sketch)
           for lines = (calculate-bond bond-node sketch)
           do (setf (lines bond-node) lines))))
-                 
+
 (defun optimize-sketch (sketch)
   ;; Set up which atoms we want to render and which ones we want to label
   (loop for atom-node in (atom-nodes sketch)
@@ -506,8 +562,18 @@ This will place the calculated bond on one or the other side of the x1,y1-x2,y2
   (loop for line in (lines node)
         do (draw-bond scene line)))
 
-
-(defvar *svg-stylesheet*
+(defmethod render-node (scene (annotation optimized-annotation))
+  (cl-svg:draw scene (:line
+                      :x1 (car (average-pos annotation)) :y1 (cdr (average-pos annotation))
+                      :x2 (car (closest-annotation-pos annotation)) :y2 (cdr (closest-annotation-pos annotation))
+                      :class "annotate"
+                      ))
+  (cl-svg:text scene (:x (car (closest-annotation-pos annotation))
+                      :y (cdr (closest-annotation-pos annotation))
+                      :class "element")
+               (text annotation)))
+
+(defparameter *svg-stylesheet*
   ".name {
      font: italic 9px sans-serif;
      fill: red;
@@ -524,6 +590,12 @@ This will place the calculated bond on one or the other side of the x1,y1-x2,y2
      stroke-width: 2;
      stroke-linecap: round;
    }
+   .annotate {
+     fill: none;
+     stroke: blue;
+     stroke-width: 1;
+     stroke-linecap: round;
+   }
    .dash {
      stroke-dasharray: 5, 5;
    }
@@ -566,7 +638,10 @@ This will place the calculated bond on one or the other side of the x1,y1-x2,y2
         do (render-dbg "bond-node #~a~%" count)
         do (render-node scene bond-node))
   (loop for atom-node in (atom-nodes sketch)
-        do (render-node scene atom-node)))
+        do (render-node scene atom-node))
+  (loop for annotation in (optimized-annotations sketch)
+        do (render-node scene annotation)))
+
 
 
 (defun svg (sketch2d &key (toplevel t) (width 1000) (xbuffer 0.1) (ybuffer 0.1) before-render after-render (id "") show-names (scale 5) (margin 40))
@@ -588,12 +663,21 @@ The caller provided functions should use cl-svg to render additional graphics."
                (geom:vec x y 0.0))))
       (let ((sketch (generate-sketch sketch2d 0 0)))
         (layout-sketch sketch #'transform-point)
+        (reinitialize-instance sketch
+                               :optimized-annotations (optimize-annotations sketch2d (atom-nodes sketch) (atoms-to-nodes sketch)))
+        (loop for optimized-annotation in (optimized-annotations sketch)
+              for pos = (closest-annotation-pos optimized-annotation)
+              do (setf
+                  xmin (min xmin (car pos))
+                  ymin (min ymin (cdr pos))
+                  xmax (max xmax (car pos))
+                  ymax (max ymax (cdr pos))))
         (setf xmin (floor (- xmin margin))
               ymin (floor (- ymin margin))
               xmax (ceiling (+ xmax margin))
               ymax (ceiling (+ ymax margin))
-              (width sketch) (- xmax xmin)
-              (height sketch) (- ymax ymin))
+              (width sketch) (* 1.2 (- xmax xmin))
+              (height sketch) (* 1.2 (- ymax ymin)))
         (optimize-sketch sketch)
         (render-dbg "About to render-sketch  bond-nodes ~a~%" (length (bond-nodes sketch)))
         ;; It looks like for jupyter we need to specify width and height

src/lisp/sketch2d/sketch2d.lisp

@@ -16,6 +16,9 @@ I'll use SMARTS to identify features.
 I'll use an angle term instead of a bond term.
 |#
 
+(defclass annotation ()
+  ((text :initarg :text :reader text)
+   (atom-names :initarg :atom-names :reader atom-names)))
 
 (defclass sketch2d ()
   ((molecule :initarg :molecule :accessor molecule)
@@ -27,7 +30,8 @@ I'll use an angle term instead of a bond term.
    (chiral-infos :initform nil :initarg :chiral-infos :accessor chiral-infos)
    (dynamics :initarg :dynamics :accessor dynamics)
    (sketch-atoms-to-original :initarg :sketch-atoms-to-original :accessor sketch-atoms-to-original)
-   (original-molecule :initarg :original-molecule :accessor original-molecule)))
+   (original-molecule :initarg :original-molecule :accessor original-molecule)
+   (annotations :initform nil :initarg :annotations :reader annotations)))
 
 (defparameter *max-tries* 3)
 (defconstant +oozp-scale+ 1.0)
@@ -1229,8 +1233,7 @@ to check if two line segments (bonds) overlap/intersect
 
 (defun sketch2d-molecule (molecule &key accumulate-coordinates verbose)
   (let* ((sketch (setup-simulation molecule :accumulate-coordinates accumulate-coordinates :verbose verbose))
-         (dynamics (dynamics sketch))
-         #+(or)(scoring-function (dynamics:scoring-function dynamics)))
+         (dynamics (dynamics sketch)))
     (sketch2d-dynamics sketch dynamics :accumulate-coordinates accumulate-coordinates)
     sketch))
 
@@ -1270,15 +1273,15 @@ are in the order (low, middle, high) and the column eigen-vectors are in the sam
         (chem:apply-transform-to-atoms molecule transposed-transform)))))
 
 
-(defgeneric do-sketch2d (thing &key accumulate-coordinates max-tries verbose one-shot))
+(defgeneric do-sketch2d (thing &key accumulate-coordinates max-tries verbose fixup))
 
-(defmethod do-sketch2d (molecule &key accumulate-coordinates (max-tries *max-tries*) (one-shot nil) verbose)
+(defmethod do-sketch2d (molecule &key accumulate-coordinates (max-tries *max-tries*) verbose fixup )
   (let ((count 1)
         problems
         unfrozen
         sketch)
     (setf sketch (sketch2d-molecule molecule :accumulate-coordinates accumulate-coordinates :verbose verbose))
-    (when one-shot (return-from do-sketch2d sketch))
+    (unless fixup (return-from do-sketch2d sketch))
     (loop named problem-loop
           for try-index below max-tries
           do (multiple-value-setq (problems unfrozen)
@@ -1310,10 +1313,10 @@ are in the order (low, middle, high) and the column eigen-vectors are in the sam
     (insert (car lst) (insertion-sort (cdr lst) compare key) compare key)))
 
 
-(defmethod do-sketch2d ((aggregate chem:aggregate) &key accumulate-coordinates (max-tries *max-tries*) verbose one-shot)
+(defmethod do-sketch2d ((aggregate chem:aggregate) &key accumulate-coordinates (max-tries *max-tries*) verbose fixup)
   (if (= (chem:content-size aggregate) 1)
       (do-sketch2d (chem:content-at aggregate 0) :accumulate-coordinates accumulate-coordinates :max-tries max-tries :verbose verbose
-        :one-shot one-shot)
+        :fixup fixup)
       (error "sketch2d only accepts a molecule or an aggregate with a single molecule")))
 
 #+debug-sketch2d
@@ -1500,12 +1503,35 @@ are in the order (low, middle, high) and the column eigen-vectors are in the sam
                  (t (warn "Handle config ~a for atom ~a" config chiral-atom)))))))
 
 
+(defgeneric sketch2d (molecule &key accumulate-coordinates verbose fixup show-names&allow-other-keys))
 
-(defun sketch2d (molecule &key accumulate-coordinates verbose one-shot use-structure show-names)
+(defmethod sketch2d ((obj topology:topology) &key accumulate-coordinates verbose fixup show-names &allow-other-keys)
+  (let* ((*show-names* (or *show-names* show-names))
+         (molecule (topology:build-one-molecule-for-topology obj))
+         (sketch2d (do-sketch2d molecule :accumulate-coordinates accumulate-coordinates :verbose verbose :fixup fixup))
+         )
+    (multiple-value-bind (stereochemistry-types configurations cips)
+        (chem:calculate-stereochemistry molecule :use-structure t)
+      (augment-sketch-with-stereochemistry nil sketch2d cips stereochemistry-types configurations)
+      (let* ((plugs (topology:plugs obj))
+             (annotations (loop for plug in (alexandria:hash-table-values plugs)
+                                for plug-name = (topology:name plug)
+                                for plug-bonds = (topology:plug-bonds plug)
+                                for plug-atom-names = (loop for plug-bond across plug-bonds
+                                                            collect (topology:atom-name plug-bond))
+                                collect (make-instance 'annotation
+                                                       :text (format nil "~s" plug-name)
+                                                       :atom-names plug-atom-names))))
+        (reinitialize-instance sketch2d :annotations annotations)
+    sketch2d))))
+
+
+(defmethod sketch2d ((molecule chem:molecule) &key accumulate-coordinates verbose fixup use-structure show-names &allow-other-keys)
   "Calculate a sketch2d from the molecule.  If use-structure is set then calculate the stereochemical configuration of atoms
 using the positions of the atoms in the molecule, otherwise get the configuration from the _Configuration slot."
+  (declare (optimize (debug 3)))
   (let ((*show-names* (or *show-names* show-names))
-        (sketch2d (do-sketch2d molecule :accumulate-coordinates accumulate-coordinates :verbose verbose :one-shot one-shot)))
+        (sketch2d (do-sketch2d molecule :accumulate-coordinates accumulate-coordinates :verbose verbose :fixup fixup)))
     (multiple-value-bind (stereochemistry-types configurations cips)
         (chem:calculate-stereochemistry molecule :use-structure use-structure)
       (augment-sketch-with-stereochemistry use-structure sketch2d cips stereochemistry-types configurations)
@@ -1578,7 +1604,7 @@ Otherwise pass a function that takes two atoms and returns T if they are matchab
               (chem:load-coordinates-into-vector energy-function (dynamics:coordinates dynamics))
               (sketch2d-dynamics new-sketch dynamics :unfrozen unfrozen :unfreeze nil
                                                      :accumulate-coordinates accumulate-coordinates)
-              (augment-sketch-with-stereochemistry use-structure new-sketch cips sterechemistry-types configurations)
+              (augment-sketch-with-stereochemistry use-structure new-sketch cips stereochemistry-types configurations)
               #+debug-sketch2d(setf (debug-info new-sketch) debug-info)
               new-sketch)))))))
 

src/lisp/topology/jupyter.lisp

@@ -1,10 +1,5 @@
 (in-package :topology-jupyter)
 
-(defun sketch-svg (topology &rest args)
-  "Generate an svg sketch of the topology - send args to sketch2d"
-  (let ((mol (topology:build-one-molecule-for-topology topology)))
-    (sketch2d:svg (apply 'sketch2d:sketch2d mol args))))
-
 (defclass node ()
   ((label :initarg :label :reader label)
    (uid :initarg :uid :reader uid)
@@ -152,11 +147,11 @@
     (build-cytoscape elements)))
 
 (defmethod cando-widgets::show-on-pane (pane-instance (object topology:topology) &rest rest &key &allow-other-keys)
-  (apply 'cando-widgets::show-on-pane pane-instance (topology-jupyter:sketch-svg object) rest))
+  (apply 'cando-widgets::show-on-pane pane-instance (sketch2d:svg (sketch2d:sketch2d object)) rest))
 
 (defmethod cando-widgets::show-on-pane (pane-instance (object symbol) &rest rest &key &allow-other-keys)
   (let ((top (chem:find-topology object)))
-    (apply 'cando-widgets::show-on-pane pane-instance (topology-jupyter:sketch-svg top) rest)))
+    (apply 'cando-widgets::show-on-pane pane-instance top rest)))
 
 
 (defmethod cando-widgets::show-on-pane (pane-instance (object topology:oligomer-shape) &rest rest &key &allow-other-keys)

src/lisp/topology/shape.lisp

@@ -562,7 +562,7 @@ callback-sidechain-rotamer-indexes - A lambda that takes the lambda-list (oligom
 
 (defun build-all-oligomer-shapes-from-internals-in-coordinates (assembler coords)
   (loop for oligomer-shape in (topology:oligomer-shapes assembler)
-        do (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape)))
+        do (build-atom-tree-external-coordinates-and-adjust assembler coords oligomer-shape oligomer-shape)))
 
 (defun random-oligomer-shape-aggregate (oligomer-shape)
   "Generate a random oligomer-shape and return the aggregate"

src/lisp/topology/steppers.lisp

@@ -598,9 +598,8 @@ This is an vector of rotamer-index values enumerated from 0...(number-of-rotamer
   (let ((rotamer-state (rotamer-state permissible-rotamers index)))
       (write-rotamers permissible-rotamers rotamer-state)))
 
-(defun first-rotamers (oligomer-shape permissible-rotamers)
+(defun first-rotamers (permissible-rotamers)
   "Return a vector of the first rotamer index for each position"
-  (declare (ignore oligomer-shape))
   (let ((len (length (permissible-rotamer-vector permissible-rotamers))))
     (make-instance 'rotamer-indexes
                    :rotamer-indexes (make-array len :initial-contents
@@ -610,9 +609,8 @@ This is an vector of rotamer-index values enumerated from 0...(number-of-rotamer
                                                       :for index = (aref (allowed-rotamer-indexes permissible-rotamer) rnd)
                                                       :collect index)))))
 
-(defun random-rotamers (oligomer-shape permissible-rotamers)
+(defun random-rotamers (permissible-rotamers)
   "Return a vector of random rotamer indexes"
-  (declare (ignore oligomer-shape))
   (let ((len (length (permissible-rotamer-vector permissible-rotamers))))
     (make-instance 'rotamer-indexes
                    :rotamer-indexes (make-array len :initial-contents