[gambit-list] JavaScript backend

Bradley Lucier lucier at math.purdue.edu
Tue Jan 12 23:39:08 EST 2016


> On Jan 12, 2016, at 11:11 PM, Marc Feeley <feeley at iro.umontreal.ca> wrote:
> 
> What is the transformation you have in mind?  It seems to me that for the optimization to operate at the basic block level, the construction of the vector and the indexing would have to be very close to each other, rather unlikely in complex code.

I think of transforming angle:apply-force in:

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; safe versions of some routines

(define (safe-sqrt x)
  (sqrt (max x 0.0)))

(define (safe-acos x)
  (acos (max -1.0 (min 1.0 x))))

; Euclidean (three-space) vector constructor, accessors, and modifiers

(define (make-euclidean x y z)
  (f64vector x y z))

(define (euclidean-x v) (f64vector-ref v 0))
(define (euclidean-y v) (f64vector-ref v 1))
(define (euclidean-z v) (f64vector-ref v 2))

(define (set-euclidean-x! v a) (f64vector-set! v 0 a))
(define (set-euclidean-y! v a) (f64vector-set! v 1 a))
(define (set-euclidean-z! v a) (f64vector-set! v 2 a))

;(define euclidean-zero '#f64(0.0 0.0 0.0))

; structure of atoms, bonds, angles, and torsions
; generally speaking, properties includes (shared) values
; specific to each atom (bond, angle, ...) type

(define (make-atom-properties mass charge)
  (f64vector mass charge))

(define (atom-properties-charge properties)
  (f64vector-ref properties 1))

(define (atom-properties-mass properties)
  (f64vector-ref properties 0))

(define (make-atom position velocity force index properties)
  (vector position velocity force index properties))

(define (atom-position atom)
  (vector-ref atom 0))

(define (atom-velocity atom)
  (vector-ref atom 1))

(define (atom-force atom)
  (vector-ref atom 2))

(define (atom-index atom)
  (vector-ref atom 3))

(define (atom-properties atom)
  (vector-ref atom 4))

(define (atom-mass atom)
  (atom-properties-mass (atom-properties atom)))

(define (atom-charge atom)
  (atom-properties-charge (atom-properties atom)))

(define (make-angle atom1 atom2 atom3 ideal-angle force-constant)
  (vector atom1 atom2 atom3 ideal-angle force-constant))

(define (angle-atom1 angle)
  (vector-ref angle 0))

(define (angle-atom2 angle)
  (vector-ref angle 1))

(define (angle-atom3 angle)
  (vector-ref angle 2))

(define (angle-ideal-angle angle)
  (vector-ref angle 3))

(define (angle-force-constant angle)
  (vector-ref angle 4))

(define (make-bond atom1 atom2 ideal-length force-constant)
  (vector atom1 atom2 ideal-length force-constant))

(define (bond-atom1 bond)
  (vector-ref bond 0))

(define (bond-atom2 bond)
  (vector-ref bond 1))

(define (bond-ideal-length bond)
  (vector-ref bond 2))

(define (bond-force-constant bond)
  (vector-ref bond 3))

; Euclidean vector operations

(define (euclidean:copy! v1 v2)
  (set-euclidean-x! v1 (euclidean-x v2))
  (set-euclidean-y! v1 (euclidean-y v2))
  (set-euclidean-z! v1 (euclidean-z v2)))

;(define (euclidean:= v1 v2)
;  (and (= (euclidean-x v1) (euclidean-x v2))
;       (= (euclidean-y v1) (euclidean-y v2))
;       (= (euclidean-z v1) (euclidean-z v2))))

;(define (euclidean:zero? v)
;  (euclidean:= v euclidean-zero))

;(define (euclidean:zero! v)
;  (euclidean:copy! v euclidean-zero))

(define (euclidean:add v1 v2)
  (make-euclidean (+ (euclidean-x v1) (euclidean-x v2))
		  (+ (euclidean-y v1) (euclidean-y v2))
		  (+ (euclidean-z v1) (euclidean-z v2))))


(define (euclidean:subtract v1 v2)
  (make-euclidean (- (euclidean-x v1) (euclidean-x v2))
		  (- (euclidean-y v1) (euclidean-y v2))
		  (- (euclidean-z v1) (euclidean-z v2))))

(define (euclidean:dot v1 v2)
  (+ (* (euclidean-x v1) (euclidean-x v2))
     (* (euclidean-y v1) (euclidean-y v2))
     (* (euclidean-z v1) (euclidean-z v2))))

(define (euclidean:cross v1 v2)
  (make-euclidean (- (* (euclidean-y v1) (euclidean-z v2)) 
		     (* (euclidean-z v1) (euclidean-y v2)))
		  (- (* (euclidean-z v1) (euclidean-x v2))
		     (* (euclidean-x v1) (euclidean-z v2)))
		  (- (* (euclidean-x v1) (euclidean-y v2))
		     (* (euclidean-y v1) (euclidean-x v2)))))

(define (euclidean:length v)
  (sqrt (euclidean:dot v v)))

(define (euclidean:scale a v)
  (make-euclidean (* a (euclidean-x v))
		  (* a (euclidean-y v))
		  (* a (euclidean-z v))))

(define (euclidean:perpendicular-component v1 v2)
  (euclidean:subtract v1 
		      (euclidean:scale (/ (euclidean:dot v1 v2) 
					  (euclidean:dot v2 v2))
				       v2)))

(define (euclidean:angle v1 v2)
  (safe-acos (/ (euclidean:dot v1 v2) 
	   (* (euclidean:length v1) (euclidean:length v2)))))


(define (angle:derivative-of-energy angle theta)
  (let ((tdif (- theta (angle-ideal-angle angle))))
    (* (angle-force-constant angle) (* tdif (+ 1.0 (* 1.508 tdif tdif))))))


(define (angle:apply-force angle)
  (let* ((atom1 (angle-atom1 angle))
	 (atom2 (angle-atom2 angle))
	 (atom3 (angle-atom3 angle))
	 (r1r2 (euclidean:subtract (atom-position atom1) (atom-position atom2)))
	 (r3r2 (euclidean:subtract (atom-position atom3) (atom-position atom2)))
	 (length-r1r2^2 (euclidean:dot r1r2 r1r2))
	 (length-r3r2^2 (euclidean:dot r3r2 r3r2))
	 (length-r1r2 (sqrt length-r1r2^2))
	 (length-r3r2 (sqrt length-r3r2^2))
	 (cos-theta (/ (euclidean:dot r1r2 r3r2)
		       (* length-r1r2 length-r3r2)))
	 (theta (safe-acos cos-theta))
	 (derivative-of-energy (/ (* 2.0 
				     (angle-force-constant angle) 
				     (- theta (angle-ideal-angle angle)))
				  (sin theta)))
	 (F1 (euclidean:subtract (euclidean:scale (/ (* cos-theta derivative-of-energy)
						     length-r1r2^2)
						  r1r2)
				 (euclidean:scale (/ derivative-of-energy
						     (* length-r1r2 length-r3r2))
						  r3r2)))
	 (F3 (euclidean:subtract (euclidean:scale (/ (* cos-theta derivative-of-energy)
						     length-r3r2^2)
						  r3r2)
				 (euclidean:scale (/ derivative-of-energy
						     (* length-r1r2 length-r3r2))
						  r1r2))))
    (euclidean:copy! (atom-force atom1) (euclidean:add (atom-force atom1) F1))
    (euclidean:copy! (atom-force atom2) (euclidean:subtract
					 (euclidean:subtract (atom-force atom2) 
							     F1)
					 F3))
    (euclidean:copy! (atom-force atom3) (euclidean:add (atom-force atom3) F3))))

(define (bond:compute-length-force bond)
  (let ((atom1 (bond-atom1 bond))
	(atom2 (bond-atom2 bond)))
    (let ((r1r2 (euclidean:subtract (atom-position atom1) 
				    (atom-position atom2))))
      (let ((length-of-r1r2 (euclidean:length r1r2)))
	(let ((derivative-of-energy (* (bond-force-constant bond) 
				       (- length-of-r1r2 (bond-ideal-length bond)))))
	  (euclidean:scale (/ derivative-of-energy
			      length-of-r1r2)
			   r1r2))))))

(define (bond:apply-length-force bond)
  (let ((atom1 (bond-atom1 bond))
	(atom2 (bond-atom2 bond))
	(force (bond:compute-length-force bond)))
    (euclidean:copy! (atom-force atom1) 
		     (euclidean:subtract (atom-force atom1)
					 force))
    (euclidean:copy! (atom-force atom2) 
		     (euclidean:add (atom-force atom2)
				    force))))

(define (non-bond:compute-vdw-force atom1 atom2 one-four)
  (let* ((r1r2 (euclidean:subtract (atom-position atom1) 
				   (atom-position atom2)))
	 (length-of-r1r2 (euclidean:length r1r2))
	 (vdw-radius (if one-four
			 (+ (atom-half-vdw-1-4-radius atom1) 
			    (atom-half-vdw-1-4-radius atom2))
			 (+ (atom-half-vdw-radius atom1) 
			    (atom-half-vdw-radius atom2))))
	 (vdw-const (if one-four
			(* (atom-sqrt-vdw-1-4-const atom1)
			   (atom-sqrt-vdw-1-4-const atom2))
			(* (atom-sqrt-vdw-const atom1)
			   (atom-sqrt-vdw-const atom2))))
	 (ratio (/ vdw-radius
		   length-of-r1r2))
	 (ratio^2 (* ratio ratio))
	 (ratio^6 (* ratio^2 ratio^2 ratio^2))
	 (derivative-of-energy (* vdw-const
				  ratio
				  ratio^6
				  (- 1.0 (* 2.0 ratio^6)))))
    (euclidean:scale (/ derivative-of-energy
			length-of-r1r2)
		     r1r2)))

(define (non-bond:apply-vdw-force atom1 atom2 one-four)
  (let ((force (non-bond:compute-vdw-force atom1 atom2 one-four)))
    (euclidean:copy! (atom-force atom1) 
		     (euclidean:subtract (atom-force atom1)
					 force))
    (euclidean:copy! (atom-force atom2) 
		     (euclidean:add (atom-force atom2)
				    force))))

(define (non-bond:compute-electrostatic-force atom1 atom2 one-four)
  (let* ((r1r2 (euclidean:subtract (atom-position atom1) 
				   (atom-position atom2)))
	 (length-of-r1r2^2 (euclidean:dot r1r2 r1r2))
	 (length-of-r1r2 (sqrt length-of-r1r2^2))
	 (epsilon (if one-four
		      non-bond:electrostatic-1-4-scaling
		      1.0)))
    (euclidean:scale (/ (* non-bond:electrostatic-constant 
			   epsilon
			   (atom-charge atom1) 
			   (atom-charge atom2))
			(* non-bond:dielectric-constant length-of-r1r2^2 length-of-r1r2))
		     r1r2)))

(define (non-bond:apply-electrostatic-force atom1 atom2 one-four)
  (let ((force (non-bond:compute-electrostatic-force atom1 atom2 one-four)))
    (euclidean:copy! (atom-force atom1) 
		     (euclidean:subtract (atom-force atom1)
					 force))
    (euclidean:copy! (atom-force atom2) 
		     (euclidean:add (atom-force atom2)
				    force))))

(define (non-bond:apply-vdw-and-electrostatic-force atom1 atom2 one-four)
  (let ((force (euclidean:add (non-bond:compute-electrostatic-force atom1 atom2 one-four)
			      (non-bond:compute-vdw-force atom1 atom2 one-four))))
    (euclidean:copy! (atom-force atom1) 
		     (euclidean:subtract (atom-force atom1)
					 force))
    (euclidean:copy! (atom-force atom2) 
		     (euclidean:add (atom-force atom2)
				    force))))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

into

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

(define angle:apply-force
  (lambda (angle)
    (let* ((g287 (vector-ref angle 0))
           (g288 (vector-ref g287 2))
           (g289 (f64vector-ref g288 0))
           (g290 (vector-ref g287 0))
           (g291 (f64vector-ref g290 0))
           (g292 (vector-ref angle 1))
           (g293 (vector-ref g292 0))
           (g294 (f64vector-ref g293 0))
           (g295 (- g291 g294))
           (g296 (vector-ref angle 2))
           (g297 (vector-ref g296 0))
           (g298 (f64vector-ref g297 0))
           (g299 (- g298 g294))
           (g300 (* g295 g299))
           (g301 (f64vector-ref g290 1))
           (g302 (f64vector-ref g293 1))
           (g303 (- g301 g302))
           (g304 (f64vector-ref g297 1))
           (g305 (- g304 g302))
           (g306 (* g303 g305))
           (g307 (f64vector-ref g290 2))
           (g308 (f64vector-ref g293 2))
           (g309 (- g307 g308))
           (g310 (f64vector-ref g297 2))
           (g311 (- g310 g308))
           (g312 (* g309 g311))
           (g313 (+ g300 g306 g312))
           (g314 (* g295 g295))
           (g315 (* g303 g303))
           (g316 (* g309 g309))
           (g317 (+ g314 g315 g316))
           (g318 (sqrt g317))
           (g319 (* g299 g299))
           (g320 (* g305 g305))
           (g321 (* g311 g311))
           (g322 (+ g319 g320 g321))
           (g323 (sqrt g322))
           (g324 (* g318 g323))
           (g325 (/ g313 g324))
           (g326 (vector-ref angle 4))
           (g327 (min 1. g325))
           (g328 (max -1. g327))
           (g329 (acos g328))
           (g330 (vector-ref angle 3))
           (g331 (- g329 g330))
           (g332 (* 2. g326 g331))
           (g333 (sin g329))
           (g334 (/ g332 g333))
           (g335 (* g325 g334))
           (g336 (/ g335 g317))
           (g337 (* g336 g295))
           (g338 (/ g334 g324))
           (g339 (* g338 g299))
           (g340 (- g337 g339))
           (g341 (+ g289 g340))
           (g342 (f64vector-set! g288 0 g341))
           (g343 (f64vector-ref g288 1))
           (g344 (* g336 g303))
           (g345 (* g338 g305))
           (g346 (- g344 g345))
           (g347 (+ g343 g346))
           (g348 (f64vector-set! g288 1 g347))
           (g349 (f64vector-ref g288 2))
           (g350 (* g336 g309))
           (g351 (* g338 g311))
           (g352 (- g350 g351))
           (g353 (+ g349 g352))
           (g354 (f64vector-set! g288 2 g353))
           (g355 (vector-ref g292 2))
           (g356 (f64vector-ref g355 0))
           (g357 (- g356 g340))
           (g358 (/ g335 g322))
           (g359 (* g358 g299))
           (g360 (* g338 g295))
           (g361 (- g359 g360))
           (g362 (- g357 g361))
           (g363 (f64vector-set! g355 0 g362))
           (g364 (f64vector-ref g355 1))
           (g365 (- g364 g346))
           (g366 (* g358 g305))
           (g367 (* g338 g303))
           (g368 (- g366 g367))
           (g369 (- g365 g368))
           (g370 (f64vector-set! g355 1 g369))
           (g371 (f64vector-ref g355 2))
           (g372 (- g371 g352))
           (g373 (* g358 g311))
           (g374 (* g338 g309))
           (g375 (- g373 g374))
           (g376 (- g372 g375))
           (g377 (f64vector-set! g355 2 g376))
           (g378 (vector-ref g296 2))
           (g379 (f64vector-ref g378 0))
           (g380 (+ g379 g361))
           (g381 (f64vector-set! g378 0 g380))
           (g382 (f64vector-ref g378 1))
           (g383 (+ g382 g368))
           (g384 (f64vector-set! g378 1 g383))
           (g385 (f64vector-ref g378 2))
           (g386 (+ g385 g375)))
      (f64vector-set! g378 2 g386))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

etc.  After inlining all the basic operations are close to one another.

Brad


More information about the Gambit-list mailing list