Last active
December 14, 2015 22:48
-
-
Save jsmpereira/5160610 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
(defun build-matrix () | |
(let* ((n-nodes (length (nodes))) | |
(dim `(,n-nodes ,n-nodes))) | |
(gsll:set-zero (grid:make-foreign-array 'double-float :dimensions dim)))) | |
(defparameter *adj-matrix* nil) | |
(defparameter *diag-matrix* nil) | |
(defmacro with-matrix (matrix &body body) | |
`(loop for i below (first (grid:dimensions ,matrix)) do | |
(loop for j below (nth 1 (grid:dimensions ,matrix)) | |
,@body))) | |
(defun adj-matrix (adj-matrix) | |
(with-matrix adj-matrix do | |
(let ((node (node-by-id i))) | |
(loop for neighbour across (neighbours node) | |
when (= (address neighbour) j) do | |
(setf (grid:aref adj-matrix j i) 1.0d0))))) | |
(defun vector-sum (adj-matrix) | |
(let ((sum (grid:make-foreign-array 'double-float :dimensions (first (grid:dimensions adj-matrix))))) | |
(with-matrix adj-matrix do | |
(setf (grid:aref sum i) (incf (grid:aref sum i) (grid:aref adj-matrix j i)))) | |
sum)) | |
(defun diag-from-vector-sum (diag-matrix vector-sum) | |
(with-matrix diag-matrix do | |
(when (= i j) | |
(setf (grid:aref diag-matrix i j) (grid:aref vector-sum i))))) | |
(defun laplacian-matrix (adj-matrix diag-matrix) | |
(gsll:elt- diag-matrix adj-matrix)) | |
(defun fiedler-value () | |
"We want the second eigen value." | |
(setf *adj-matrix* (build-matrix)) | |
(setf *diag-matrix* (build-matrix)) | |
(adj-matrix *adj-matrix*) | |
(diag-from-vector-sum *diag-matrix* (vector-sum *adj-matrix*)) | |
(gsll:eigenvalues (laplacian-matrix *adj-matrix* *diag-matrix*))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
A sane return value (I'd wager) for (gsll:eigenvalues):
Sometimes I get stuff like this: