hamf.lisp

Example of using GSL inside lisp

(clc:clc-require :gsll)
(clc:clc-require :cgn)
;(load "/home/faguayo/INSTALLS/MATLISP/matlisp-2_0beta-2003-10-14/start.lisp")

(defpackage :myown
 (:use :common-lisp
	:gsll
	:cgn))

;(defpackage :myown
; (:use :common-lisp
;       :matlisp
;       :cgn))

(in-package :myown)

;;;;;;;;;;;;;;;;;;;;; Funciones ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;toma una lista '((x y) (x y)...) y la grafica
;; si se especifica un filename, deja guardado el archivo de datos
(defun gnuplotea-2d ( lista &key filename )
  (let ( (file (string ".cgn.dat")) )
    (if filename (setf file filename) )
    (with-open-file (stream file :direction :output  :if-does-not-exist :create :if-exists :overwrite )
      (dotimes (i (list-length lista))
	(format stream "~D   ~D   ~%" (car (nth i lista)) (cadr (nth i lista)) )
	)
      )
  (with-gnuplot ( 'linux )
  (format-gnuplot "plot ~s u 1:2 w l" file)
  (print-graphic))
  )
  )


;;toma un 'marray' de dimension 2 y hace un grafico de la superficie (los primeros nnx y nny puntos)
;; sin interpolar
;; si se especifica un filename, deja guardado el archivo de datos
(defun gnuplotea-3d ( matriz nnx nny &key filename )
  (let ( (file (string ".cgn.dat")) )
    (if filename (setf file filename) )
    (with-open-file (stream file :direction :output  :if-does-not-exist :create :if-exists :supersede )
      (dotimes (i nnx)
	(dotimes (j nny)
	  (format stream "~D   ~D   ~D   ~%" i j (maref matriz i j) )
	  )
	(format stream "~%")
	)
      )
    (with-gnuplot ( 'linux )
      (format-gnuplot "set pm3d corners2color c1 map")
      (format-gnuplot "set size ratio -1")
      (format-gnuplot "splot ~s u 1:2:3 w pm3d" file)
      (print-graphic))
    ) 
  )

;; escribe una matriz con formato humano
;(defun escribe (matriz)
;  (dotimes (yi *Ny*)
;       (dotimes (xi *Nx*)
;	    (format t " ~6$ " (maref matriz xi yi))
;	    )
;       (format t "~%")
;       )
;  )

;;;;;;;;;;;;;;;;;;; FIN FUNCIONES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;; OOP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; definimos la clase "sitioa" por sitio activo
;; tiene 2 slots
;; - coordenadas (i,j)
;; - una lista con sus vecinos, que son instancias de la misma clase
;; mas adelante podriamos definirle una "cantidad de combustible restante" y una
;; "radiacion recibida" o temperatura.
(defclass sitioa ()
  (
   (posicion
    :initarg :posicion
    :reader posicion)
   (vecinos
    :initform '() 
    :reader vecinos)
   (orden
    :initarg :orden
    :reader orden)
   )
  )

;; definimos un metodo/funcion que dada una lista y un sitio de esta, busca que elementos dentro de la misma lista 
;; son los "l" (lx,ly definen la elipse) vecinos de este elemento.
;; como es una relacion "biyectiva" podriamos definirlos a ambos como vecinos
;; simultaneamente, pero aun no se como evitar que la lista se repita.

;; primero un metodo que nos diga si un sitio (sitio2) esta dentro de la zona
;; de influencia del sitio1
(defgeneric en-zona (sitio1 sitio2 lx ly)
  (:documentation "Ve si sitio2 esta dentro de la zona de influencia del sitio1. 
lx y ly son los parametros de la elipse"))

(defmethod en-zona ((sitio1 sitioa) (sitio2 sitioa) lx ly)
  (let* (
	 (v (mapcar #'- (posicion sitio1) (posicion sitio2)) )
	 (a (car  v))
	 (b (cadr v))
	 )
    (if (<= (+ (expt (/ a lx) 2) (expt (/ b ly) 2)) 1.001) t nil)
    )
  )

;; ahora un metodo que al sitio1 le asigna como vecino el sitio2 
(defgeneric nuevo-vecino (local visita)
  (:documentation "asigna al sitioa local el sitioa visita como vecino"))
 
(defmethod nuevo-vecino ((local sitioa) (visita sitioa))
  (let* ((vec-local (slot-value local 'vecinos))
	 (ya-esta nil))
    (dotimes (i (list-length vec-local))       ;revisamos si ya esta el sitio dentro de los vecinos
      (if (equal (posicion (nth i vec-local)) (posicion visita))
	  (setf ya-esta t)) 
      )
  (unless (or ya-esta (equal (posicion local) (posicion visita)) ) ;si no esta, y si el vecino no tiene la misma posicion que el local, lo agregamos
    (if (null vec-local)                                           ;agregamos en dos casos: si la lista esta vacia la creamos, de otra forma concatenamos
      (push visita (slot-value local 'vecinos))                    ;asi tine sentido hablar del vecino del vecino del vecino... 
      (nconc (slot-value local 'vecinos) (list visita))
      )
    )
  ))
 

;;;;;;;;;;;;;;;;;;; FIN OOP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;; Programa Principal ;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defparameter *Nx* 14)       ;numero de puntos en X
(defparameter *Ny* 14)       ;numero de puntos en y
(defparameter *Pb* 0.5)     ;probabilidad del dopaje
(defvar *rng*)              ;generador de numeros aleatorios
(defvar *Arnd*)             ;Matriz con puntos randoms
(defvar *Abin*)             ;lista de 'sitiosa' que estan activos sitios activos '((x y) (x y)...)
(defvar *Aprint*)             ;como Arnd pero para imprimir
(defvar *Ham*)              ;Hamiltoniano
(defparameter J0 -1.0D0)      ;Fuerza del coupling en el hamiltoniano


(setf *rng* (make-random-number-generator +mt19937+ 0))
(setf *Arnd* (make-marray 'double-float :dimensions (list *Nx* *Ny*)))
(setf *Aprint* (make-marray 'double-float :dimensions (list *Nx* *Ny*)))
(setf *Abin* '())

(defvar *Arnd*)             ;Matriz con puntos randoms
(setf *Arnd* (make-marray 'double-float :dimensions (list 1500 1500)))

(dotimes (xi *Nx*) 
     (dotimes (yi *Ny*) 
       (let ( (rnd-temp (uniform *rng*)) )
	 (setf (maref *Arnd* xi yi) rnd-temp)
	 (if (< rnd-temp *Pb*) 
	     (if (null *Abin*)
		 (push (make-instance 'sitioa :posicion `(,xi ,yi) :orden (list-length *Abin*)) *Abin*) 
		 (nconc *Abin* (list (make-instance 'sitioa :posicion `(,xi ,yi) :orden (list-length *Abin*))))
)))))


(dotimes (i (list-length *Abin*))
  (dotimes (j (list-length *Abin*))
    (cond ((en-zona (nth i *Abin*) (nth j *Abin*) 1 1 )
	   (nuevo-vecino (nth i *Abin*) (nth j *Abin*))
))))

;(print *Abin*)
(defvar *Aprint1* '())
(setf *Aprint1* (make-marray 'double-float :dimensions (list *Nx* *Ny*)))
(dotimes (i (list-length *Abin*))
  (let  ( (pos (posicion (nth i *Abin*))) )
  (setf (maref *Aprint1* (car pos) (cadr pos)) 1.D0)
))
;(gnuplotea-3d *Aprint1* *Nx* *Ny*)
;(gnuplotea-2d *Abin* :filename "hola.dat")
;(gnuplotea-3d *Arnd* *Nx* *Ny* :filename "hola.dat" )

;(dotimes (i (list-length *Abin*))
;  (print (posicion (nth i *Abin*)))
;  (print (map 'list #'posicion (vecinos (nth i *Abin*))))
;  (print "---------")
;)

(defvar *Nh*)
(setf *Nh* (list-length *Abin*))
(setf *Ham* (make-marray 'double-float :dimensions (list *Nh* *Nh*)))
(dotimes (i *Nh*)
  (let  ((vec (vecinos (nth i *Abin*))) )
    (dolist (v-sel vec)
      (setf (maref *Ham* i (orden v-sel)) J0)
)))
;(gnuplotea-3d *Ham* *Nh* *Nh*)


;; FIN de la geometria, vamos a ver los autovalores/vectores
(defvar VAL)
(defvar VEC)
(defvar W)
(setf VAL (MAKE-MARRAY 'double-float :dimensions *Nh*))
(setf VEC (MAKE-MARRAY 'double-float :dimensions `(,*Nh* ,*Nh*)))
(setf W (MAKE-EIGEN-SYMMV *Nh*))
(EIGENVALUES-EIGENVECTORS *Ham* VAL VEC W)

;ordenamos los autovalores de menor a mayor
(defvar *diccionario* '())
(dotimes (i *Nh*)
  (let ((vall (maref VAL i)))
    (push `(,vall ,i) *diccionario*)
))
(setf *diccionario* (sort *diccionario* #'< :key #'car))

(print val)

(set-all *Aprint* 0.0D0)
;(dolist (elemento *diccionario*)
(dotimes (j (list-length *diccionario*)) 
  (let* ((elemento (nth j *diccionario*))
	 (ae     (car  elemento))
	 (indice (cadr elemento)))
    (cond ((< ae 0)
	   (let ((vector (column VEC indice)))
	     (dotimes (i *Nh*)
	       (let* ((vv (nth i *Abin*))
		      (vx (car  (posicion vv)))
		      (vy (cadr (posicion vv))))
		 (setf (maref *Aprint* vx vy) 
		       (+ (maref *Aprint* vx vy)
			  (expt (maref vector i) 2)))
		 )))))))

(gnuplotea-3d *Aprint* *Nx* *Ny*)
(gnuplotea-3d *Aprint1* *Nx* *Ny*)
(close-gnuplot)	

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



(format t "~1$" *Abin*)
(escribe *Arnd*)
(list-length *Abin*)







(print VAL)

(marray-dimensions a)

(setf V0 (MAKE-MARRAY 'double-float :dimensions 2))
(setf V1 (MAKE-MARRAY 'double-float :dimensions 2))
(setf V2 (MAKE-MARRAY 'double-float :dimensions 2))

(setf V0 (column VEC 0))
(setf V1 (column VEC 1))
(setf V2 (column VEC 2))

(dot V0 V1)
(dot V0 V2)
(dot V1 V2)
(dot V0 V0)
(dot V1 V1)
(dot V2 V2)

(print VAL)



(gsl-lookup "gsl_matrix_set_all")
(gsl-lookup "gsl_rng_uniform")
(gsl-lookup "gsl_blas_sdot")
(gsl-lookup "gsl_rng_alloc")
(gsl-lookup "all-random-number-generators")
(documentation #'SET-ALL 'function)
(MAKE-RANDOM-NUMBER-GENERATOR random128_glibc2)

(get-random-number "unifor")


(rng-environment-setup)
(make-random-number-generator aaa 0)
(print (uniform aaa))

(all-random-number-generators)

(documentation #'GET-VALUE 'function)
(documentation #'uniform-fixnum 'function)
(examples 'get-RANDOM-NUMBER-GENE)


(setf rng (make-random-number-generator +ranlxd1+ 5))



(save-test random-number-generators
 (let ((rng (make-random-number-generator +mt19937+ 0)))
   (loop for i from 0 to 10
         collect
         (uniform-fixnum rng 1000)))
 (let ((rng (make-random-number-generator *cmrg* 0)))
   (loop for i from 0 to 10 collect (uniform rng))))


(with-gnuplot ( 'linux )
  (format-gnuplot "plot x w l")
  (print-graphic))