Archivo de la etiqueta: julia

Lo que he aprendido: Julia, Babel y Org

Ya sabéis que me gusta llevar la contraria. También que me gusta el minimalismo y usar el mínimo de herramientas posible. Así que cuando me ha dado por explorar la programación literaria y la investigación reproducible en Julia, en lugar de aprender a utilizar los cuadernos Jupyter como una persona normal, me he puesto a jugar con mi adorado modo Org y su funcionalidad para incluir bloques de código Babel1. Supongo que es lo que tiene ser la única de tu entorno (laboral) que hace estas cosas y no tener por tanto limitaciones en la elección de herramientas. Total, como en la academia no se colabora, no se va a dar el caso de que mi elección me limite. En fin, dejo la academia, que me caliento, y os explicito los motivos:

  • No instalo más software: ya estoy programando en Julia en Emacs así que no necesito mil programas y mi ordenador durará años y años.

  • Utilizo mis atajos de teclado y configuración: como sigo usando el mismo programa, no tengo que aprender dónde están otros botoncillos. Además, tengo definido cómo exportar a pdf desde Org pasando por LaTeX en mi archivo de configuración por lo que el documento final es exactamente como yo quiero.

  • Puedo cambiar de lenguaje de programación: el sistema funciona para una pila de lenguajes de programación. Puedo incluso mezclar código de diferentes lenguajes en un único documento2.

Visto esto, parece que elegir Org es una buena idea, así que paso a contaros en qué consiste la programación literaria, cómo se aplica en Org y qué tiene de especial el caso de Julia.

La idea: escribir para humanos

La idea de la programación literaria, como tantas otras buenas ideas, fue de Donald Knuth. Consiste en cambiar de lector objetivo al escribir un programa: en lugar de una secuencia de instrucciones para una máquina salpicada por algún comentario sobre la implementación, el programa se convierte en una descripción de la lógica detrás de la implementación salpicada por algún fragmento de código. De esta manera, no es necesario descifrar qué hace el programa leyendo el código porque las decisiones de diseño y la lógica se explicitan. Se puede pensar en el programa, entonces, como en una obra literaria.

Como me dedico a los métodos de cálculo, es decir, propongo una manera de calcular algo y luego demuestro mediante un ordenador que mi manera es mejor que la manera anterior, este enfoque me interesa por dos motivos:

  • Investigación reproducible: puedo escribir un artículo científico que incluya mis datos y mi código con lo que quien lo lea puede acompañarme en el proceso, obtener mis mismos resultados y verificar si que mis conclusiones son correctas.

  • Documentación útil: puedo explicar en el propio programa mi proceso mental para implementar numéricamente un cálculo que inicialmente era analítico, incluyendo matemáticas si es preciso.

Solo nos queda responder a una pregunta: ¿cómo demonios se lee esa mezcla de texto y código sin volverse una loca? Pues mediante los procesos de weaving, que deja solo lo humano, y tangling, que deja solo lo que entiende la máquina. ¡A ver si os habíais pensado que el señor Knuth no había pensado en esto! La cuestión aquí es que como tanto la documentación o explicación y el código salen de un mismo documento, ambos crecen juntos y limitamos la típica divergencia de según la documentación esta función tiene dos variables de entrada pero según el código tiene tres.

Programación literaria en Org

Ahora que sabemos qué es lo que queremos hacer, vamos a ver cómo lo hacemos en Org. No puede ser más fácil: cuando queramos meter un cacho de código escribimos

#+BEGIN_SRC lenguaje opciones_cabecera

Código

#+END_SRC

y ya está. Impresionante. Más os voy a decir: los bloques se pueden crear con las plantillas fáciles de Org situándonos en el principio de la línea y haciendo <s TAB.

Para poder ejecutar el trocito de código necesitamos primero decirle a Org que vamos a usar el lenguaje en cuestión y que ese lenguaje sea uno de los soportados. Tan sencillo como ir al archivo de configuración y añadir elementos a la lista de lenguajes:

(org-babel-do-load-languages
(quote org-babel-load-languages)
(quote (
(LENGUAJE1 . t)
(LENGUAJE2 . t)
)))

Ahora si hacemos C-c C-c sobre el trocito, lo ejecutaremos y nos saldrá el resultado debajo. Por seguridad, preguntará si queremos ejecutar y tendremos que contestarle, si somos vagos y nos gusta ver el mundo arder podemos decirle que ejecute todo el código sin preguntar con

(setq org-confirm-babel-evaluate nil)

¡Destrucción! ¡Sí!

Un detallico sobre los resultados: los trocitos de código se ejecutan por defecto como si fueran una función (opción de cabecera :results value) y solo nos devolverán el contenido si se lo pedimos explícitamente, con un return en el caso de Python, por ejemplo. Podemos cambiar la opción a :results output y, entonces, Org nos devolverá el contenido de STDOUT. Para lenguajes interpretados, podemos combinar esta opción con :session, que abre un intérprete y envía allí el código de todos los bloques que contengan dicha opción. Es un poco lío esto, pero jugando un poco con las opciones y leyendo la docu, se entiende.

También usamos las opciones de cabecera para decidir si al crear el documento exportamos el código (:exports code), los resultados (:exports results), el código y los resultados (:exports both) o nada (:exports none); para decirle a Org qué debe hacer con los resultados (:post); o para decir si queremos solo el archivo para la máquina (:tangle ARCHIVO_DESTINO), que por defecto no nos crea.

También podemos configurar cómo exporta el código al documento final. En mi caso, como el documento final es un pdf y eso pasa por LaTeX, quiero que use listingsUTF8 para los bloques de código. Se puede configurar para minted también, claro.

;; Exportar código como listings
(require 'ox-latex)
(setq org-latex-listings t)

;; Paquetes de LaTeX para código
(setq org-latex-default-packages-alist
'((""    "listingsutf8"  t)
("usenames,dvipsnames,svgnames,table"    "xcolor"    t))))

Las cuqueces y los colorinchis no se limitan al documento final, con

(setq org-src-fontify-natively t)

también tendremos colores en nuestro Org.

Ah, por cierto, el trocito de código lo podemos modificar en un buffer especial que se abre con C-c '.

El caso de Julia

Julia es uno de los lenguajes que no tiene soporte directamente en Babel porque su autor no le dio el copyright a la FSF. Por lo tanto, aparte de añadir Julia a la lista de lenguajes que puede usar Org, necesitamos los paquetes ESS y ob-julia.

Luego, hacen falta un par de líneas extra en el archivo de configuración para decirle a ob-julia dónde está el ejecutable de Julia (yo lo tengo en el PATH y por eso no le doy la ruta entera) y decirle a Emacs que use ob-julia:

;; Código Julia en Org
(setq inferior-julia-program-name "julia") ;; nombre o ruta de ejecutable
(require 'ob-julia)

(org-babel-do-load-languages
(quote org-babel-load-languages)
(quote (
(julia . t)
)))

Ale, ya puedo jugar a programar explicándome a mí misma lo que he hecho. Ahora solo me falta aplicar estas ideas al archivo de configuración de Emacs para no romperlo nunca.

Seguiremos informando.

Referencias

Documentación de Babel

Babel: active code in Org-mode

Working with source code en el manual de Org

Introduction to Literate Programming

Emacs para ciencias del dato

Julia with Emacs Org mode

Org-mode and julia: an introduction

Literate Programming Examples

Literate programming with Org-mode

Ten simple rules for writing and sharing computational analyses in Jupyter Notebooks

Drops of Jupyter notebooks: how to keep notes in the information age


Os dejo con música:


  1. También hay un paquete específico para Julia, pero ¿me gusta a mí lo fácil? No. 
  2. Lo digo como su supiese programar en múltiples lenguajes, ¿habéis visto? 

Lo que he aprendido: Julia en Emacs

He conseguido por fin poder programar en Julia usando Emacs. Eh, que no ha sido tan fácil: algo tan sencillo como instalar el julia-mode y el ess se convirtió en un infierno. Me pasó de todo.

Primero en GNU/Linux, como no estaba usando la versión más novedosa de elementary, el Emacs de los repositorios era demasiado antiguo y no compatible con los modos que me hacían falta. Me quedaban varias opciones: (i) actualizar el sistema (lo que finalmente hice); (ii) añadir un ppa con un Emacs más moderno (bien, pero prefiero tener los paquetes de mi distro que ya he organizado alguna petarda); o (iii) compilar Emacs desde fuente (ya lo hice alguna vez y no me apetecía repetir).

Luego en Windows, donde sí tenía un Emacs lo suficientemente moderno, parece ser que Julia tiene no sé que bug y se cuelga y hay que darle a C-g para que se reviva.

En fin, que actualicé el sistema operativo y en quince minutos ya tenía un entorno para programar en Julia. Bueno, eso en mi recién estrenado Juno, en Windows me conformo con Atom, ese editor que no es un editor. Qué dura es la vida a veces.

Después de contaros mis desventuras (para eso tengo un blog), paso a resumir qué hice la vez que funcionó todo.

Julia en GNU/Linux

Julia no está en los repos, nos dejan unos binarios genéricos para que descarguemos y ejecutemos sin más. Hay que elegir unos u otros según la arquitectura de nuestro procesador. Yo, como nunca sé qué tengo lo miro así:

ondiz@slimbook:~$ lscpu
Arquitectura:                        x86_64
modo(s) de operación de las CPUs:    32-bit, 64-bit

Descargados los binarios correspondientes, es útil crear un enlace simbólico a algún lugar donde ya estemos buscando ejecutables o añadir la carpeta donde los hemos descargado al PATH, a gusto del consumidor.

Yo he elegido la primera opción, así que he hecho:

sudo ln -s RUTA_A_EJECUTABLE /usr/local/bin/julia

Cuidadín que hay que poner la ruta absoluta al ejecutable bin/julia que si no no carrula.

Julia en Emacs

Pues lo que decía al principio: para tener un entorno chachi para programar en Julia solo hace falta instalar los modos julia-mode y ess. El primero es el modo oficial para editar Julia y el segundo es un paquete que ayuda a trabajar con lenguajes de programación para estadística (ESS viene de Emacs Speaks Statistics) como R, o en este caso, Julia. Se pueden hacer otras cosas, pero esta es la más simple en mi opinión.

Una vez instalados los modos, para activar una terminal juliana inferior solo hay que hacer M-x julia. Luego ya podemos abrir un archivo en la parte superior y jugar con opciones que nos da el modo ESS.

Llevo poco con el tema y solo he memorizado un par de combinaciones útiles (miento, no las he memorizado, las escribo aquí para mirarlas en el futuro):

  • C-c C-l: carga un archivo completo, lo que sería un include("archivo").

  • C-M-x: ejecuta (me encanta esta palabra) un trozo de código en el REPL.

  • C-c C-z: cambia del script al REPL o viceversa.

Hay también una pila de comandos para gestionar errores y para acceder a la documentación que algún día controlaré. O no. También me falta echarle un ojo a imenu anywhere para que me aparezcan opciones de autocompletado en cualquier lugar. Me lo dejo de deberes.

Curiosamente, lo que más trabajo me dio fue acostumbrarme a no darle a la flecha hacia arriba para repetir el último comando en el REPL. Aquí, como las combinaciones de GNU Readline están ya pilladas, hay que usar M-p y M-n (o C-↑ y C-↓) para moverse por la historia. Comint y cosas, ya sabéis, y si no, con hacer C-h b os enteráis.

Y hala, ya tengo un entorno para programar. Contadme qué usáis vosotros para que aprendamos todos, venga.

juliaEmacs

Referencias

Página oficial de Julia

Manual del modo ess

Julia en el modo ess

imenu anywhere en GitHub


Os dejo con un grupo cuyas canciones suelo berrear en los conciertos de las fiestas de los pueblos vecinos y que tiene el mérito de que no me quedase dormida encima del libro de física de primero de carrera. Era abrirlo y bostezar, oigan.

Lo que he aprendido: scraping con Julia

Tengo la intuición de que los artículos de los periódicos son cada vez más sencillos, de que se usan frases más cortas y palabras más comunes. Es posible que se deba a que yo me culturizo con la edad y que, por lo tanto, no tengo la misma percepción de la complejidad que tenía hace unos años. Como cuando veo Saber y ganar y me sé las respuestas, bueno, las es igual exagerar un poco. O puede ser que ocurra como con Willy Fogg, que hizo que una generación entera se criara oyendo a un león en traje decir “detrás de usted, estimado caballero” o frases del estilo, mientras que los dibujos de ahora hablan más en la línea de la chavalada.

Como navego entre la fina línea que separa la cordura de la sinrazón, he pensado que podría analizar diferentes periódicos en diferentes épocas y ver si mi hipótesis (ha evolucionado de intuición a hipótesis ¿habéis visto?) es correcta. Soy una científica, qué leches. Y de ahí el título de esta entrada: cuando una quiere analizar textos, necesita primero conseguirlos y qué mejor manera que usar una técnica que no puede estar más de moda.

Antes de contar qué he hecho va un disclaimer: esto es el resultado de un par de tardes de investigación por parte de alguien que no sabe programar en general y menos en Julia en particular pero que, imbuida por el espíritu hacker, se ha puesto a darle a todos los botones hasta que ha conseguido algo que más o menos funciona. Así funciona la ingeniería, hermanos.

Lo que he entendido sobre el scrapeo

Veamos, si yo me he empanado de algo, scrapear consiste en crear un robot que navegue por una web como si fuera un humano y se descargue algún tipo de dato para luego analizarlo. Para ello, básicamente hay que localizar en el código fuente de la página el elemento en cuestión para decirle a nuestro robot a dónde debe ir y qué debe filtrar.

En mi caso, fui a la hemeroteca de El País para un día determinado, localicé la url de las noticias de ese día y, para cada una de ellas, descargué el texto. Todo ello lo hice con Julia, que se note que estoy aprendiendo.

Scrapear con Julia

Utilicé los siguiente paquetes:

  • HTTP: nos da funcionalidad cliente/servidor HTTP, lo usé para descargar el HTML de la página.
  • Gumbo: un parser HTML, convierte el HTML en un árbol formado por nodos y elementos.
  • Cascadia: un selector CSS, permite filtrar el output de Gumbo por clase o tipo de elemento, por ejemplo.
  • AbstractTrees: sirve para manejar datos de tipo árbol, útil para pasear por los árboles que genera Gumbo.

Seguramente se podría hacer con otros, pero los ejemplos que fusilé y remezclé utilizaban estos, así que para qué complicarse. Los instalé con pkg desde la terminal de Julia.

El proceso tuvo dos partes: localizar las noticias y extraer el texto de cada una de ellas.

Localizar las noticias

Para localizar las noticias, fui a la hemeroteca, elegí una fecha y me fijé en que la url a la que me mandaba era https://elpais.com/tag/fecha/AAAAMMDD donde AAAAMMDD es la fecha con los cuatro dígitos del año (AAAA), seguido de los dos dígitos del mes (MM) y los dos del día (DD).

Una vez ahí, vi que cada una de las noticias estaba dentro de un elemento cuya clase era articulo-titulo y que contenía la dirección de la noticia en cuestión. Así, si conseguía extraer la url de cada noticia, podría hacer que mi robot fuera a esa dirección y descargase el texto.

Hice lo siguiente:

date = 19851010 # Elegir una fecha
url = "https://elpais.com/tag/fecha/"*date # Montar la url
res = HTTP.get(url) # Descargar HTML
body = String(res.body) # Convertir a texto
html = parsehtml(body) # Parsear HTML
articles = eachmatch(sel".articulo-titulo", html.root); # Seleccionar los elementos que contienen las noticias

# Para cada noticia
for f in PreOrderDFS(articles) # Pasea por los elementos
   if f isa HTMLElement{:a} # Si es un enlace
       url = getattr(f,"href") # Extrae la url
       # TODO: extraer texto de cada noticia
   end
end

Faltaba extraer el texto de cada noticia.

Extraer el texto

Despues, fui a varias noticias y vi que el texto de la noticia en sí estaba dentro de un elemento de clase articulo-cuerpo. Ese elemento estaba formado por otros elementos que estaban formados por otros hasta llegar a los elementos HTMLText que contenían el texto. Solo tenía que pasearme por las hojas (los elementos del árbol que no tienen a su vez elementos) extraer su texto y empalmar los trozos, saltando una línea en el caso de que el texto viniese de un párrafo, sin saltar si no.

Escribí este código y lo introduje en el TODO anterior:

html = parsed(url) # Parsear HTML
content = eachmatch(sel".articulo-cuerpo", html.root) # Seleccionar cuerpo
texto = "" # Inicialización
for elem in Leaves(content) # Pasear por las hojas
   if elem isa HTMLText # Si el elemento es texto
   # Si el texto viene de un párrafo saltar línea
       if elem.parent isa HTMLElement{:p}
           texto = texto*text(elem)*"\n"
       else
           texto = texto*text(elem) # Unir a lo anterior
       end
    end
end

Juntando todo, extrayendo el código repetido en un par de funciones y metiendo los textos de las noticias en un array, me quedó:

using HTTP, Gumbo, Cascadia, AbstractTrees

"""
    getText(HTMLNode)
Extrae texto de objeto HTMLNode
"""
function getText(content)
texto = ""
for elem in Leaves(content)
    if elem isa HTMLText
        # Si el texto viene de un párrafo saltar línea
        if elem.parent isa HTMLElement{:p}
           texto = texto*text(elem)*"\n"
        else
           texto = texto*text(elem)
        end
    end
end
return texto
end

"""
   parsed(url)
Devuelve HTML parseado a partir de URL
"""
function parsed(url)
res = HTTP.get(url)
body = String(res.body)
html = parsehtml(body)
return html
end

"""
   scrapElPais(date)
Descarga las noticias de El País para una fecha dada como AAAAMMDD
"""
function scrapElPais(date)
url = "https://elpais.com/tag/fecha/"*date
html = parsed(url)
articles = eachmatch(sel".articulo-titulo", html.root);
# Inicializar matriz que contendrá textos
raw = String[]
for f in PreOrderDFS(articles)
    if f isa HTMLElement{:a}
       url = getattr(f,"href")
       html = parsed(url)
       content = eachmatch(sel".articulo-cuerpo", html.root)
       articleText = getText(content)
       append!(raw, [articleText])
   end
end
return raw
end
################# MAIN #######################
date = "19760504" # Fecha de las noticias
news = scrapElPais(date)

Y así es cómo se scrapea un periódico sin saber scrapear ni programar ni entender cómo funciona la web. Evidentemente, esto se puede afinar, lo iré haciendo y os contaré. Ahora voy a tratar los datos: contaré el número de palabras total, por párrafo y por frase y miraré cuántas palabras aparecen en la lista de palabras más comunes. Será divertido.

Por cierto, WordPress me revienta la indentación y no tiene resaltador de sintaxis para Julia (estoy usando la de Python), para algo más decente tenéis el snippet de Gitlab.

Referencias

Julia: Introduction to Web Scraping (PHIVOLCS’ Seismic Events)

Web scraping with Julia

Web scraping en diferentes lenguajes en Rosetta Code


Os dejo con rap en galego:

Lo que he aprendido: plano de fase con PyPlot en Julia

Todo comenzó porque estaba resolviendo numéricamente una ecuación diferencial ordinaria (EDO para los amigos) y quería ver su comportamiento para diferentes condiciones iniciales. ¿Y cuál es la mejor herramienta para eso? El plano de fase. Así que pinté uno en Matlab y luego, como yo no soy de quedarme satisfecha al terminar algo, me pregunté: ¿cómo se hará esta historieta en Julia?

Después de hurgar en la documentación, dar a todos los botones y pintar cosas sin sentido alguno, lo conseguí con PyPlot, uno de los paquetes que se pueden usar en Julia para hacer gráficos.

Os cuento cómo calculé y cómo pinté. Todo es para Julia 1.0.

Unas notas sobre PyPlot

Antes de meternos en harina, no está de más recordar que Julia permite utilizar diferentes paquetes para dibujar, están entre otros el mítico Plots; GR, que va como un tiro y el que he usado yo esta vez: PyPlot, que es una interfaz para Julia del Matplotlib de Python.

No tengo especial predilección por este paquete, simplemente encontré un ejemplo que se parecía a lo que yo quería hacer y seguí la norma básica de la ingeniería: no reinventar la rueda. Norma especialmente certera cuando eres ingeniera mecánica.

Debo añadir también que PyPlot me funcionó directamente, probablemente porque tenía Miniconda instalado y tal vez Matplotlib también. Ni idea. Estoy todo el rato probando chismes, no sé ya ni lo que hago.

En fin, los gráficos en PyPlot funcionan casi como los de Matlab excepto por un par de cosas que a mí me dejaron loca y os pongo aquí para que no os pase lo mismo:

  • Para mostrar la figura actual en una ventana necesitamos hacer gcf() (get current figure), al menos desde un script.
  • PyPlot añade por defecto las figuras encima de las que ya había. Antes había opción hold on como en Matlab pero se desaconseja su uso.
  • A veces no querremos añadir las figuras encima, para borrar anterior tenemos clf() (clear figure)

Sabiendo esto ya podemos hacer gráficos sin control.

El plano de fase

Para que todo el rollo de la ecuación diferencial sea un poco menos abstracto os voy a dar un poco de contexto: estoy estudiando un sistema de un grado de libertad, es decir, una masa atada a una pared mediante un muelle y un amortiguador.

De esta manera, el plano de fase consiste en un gráfico en el que se representa la posición de la masa en el eje horizontal y su velocidad en el vertical y se pintan unas flechillas que indican la dirección en la que saldrá el sistema al darle unas condiciones iniciales de posición y velocidad.

Nos sirve para ver si aumenta o se reduce su velocidad, tiende hacia una posición de equilibrio o no, o para hacer bonito en un paper, la finalidad última de toda nuestra investigación.

Bien, la ecuación que describe el movimiento del sistema de un grado de libertad es la siguiente:

m\ddot{x} + c\dot{x}+kx = 0

donde el puntico hace referencia a la derivada respecto al tiempo y dos punticos a la segunda derivada. Los valores m (masa), c (amortiguamiento) y k (rigidez) son características del sistema y son conocidos.

Voy a hacer una magia todo loca y convertir la ecuación de segundo orden en dos ecuaciones de primer orden. La primera va a ser simplemente el cambio de variable:

v = \dot{x}

osease: que la velocidad es la derivada de la posición respecto al tiempo.

La segunda proviene de la ecuación original y viene a decir que la segunda derivada de la posición (la aceleración) es la primera derivada de la velocidad y por lo tanto:

\dot{v} = \ddot{x}= -\frac{k}{m}x - \frac{c}{m}v

Si cogemos las dos ecuaciones, reorganizamos de tal manera que las cosas con puntico estén a la izquierda y escribimos el sistema en forma matricial tenemos lo siguiente:

\begin{bmatrix}\dot{x}\\\dot{v}\end{bmatrix} =  \begin{bmatrix}0 & 1\\-k/m& -c/m\end{bmatrix}  \begin{bmatrix}x\\v\end{bmatrix}

Todo este lío nos lleva a algo, paciencia. La cuestión es que queremos ver cómo cambia la posición y la velocidad en cada punto. ¿Y qué cacharro matemático explica cómo cambian las cosas? Exactamente, la derivada. Derivada que precisamente tenemos explícitamente representada en el sistema de ecuaciones.

¿Por qué digo eso? Si yo quiero saber cuánto valen las derivadas de mis variables en el punto (1,0) no tengo más que sustituirlo en la ecuación:

\begin{bmatrix}\dot{x}\\\dot{v}\end{bmatrix} =  \begin{bmatrix}0 & 1\\-k/m& -c/m\end{bmatrix}  \begin{bmatrix}1\\0\end{bmatrix}

como m, c y k son valores conocidos, puedo obtener los valores de las derivadas. Sabiendo esto, puedo pintar una flechillas en el plano que indiquen hacia dónde tiende el sistema si le doy determinadas condiciones iniciales.

Antes de pintar las flechillas, os dejo aquí la función que he escrito en Julia para calcular las derivadas de la posición y el desplazamiento en diferentes puntos del plano. Como es una función genérica, he utilizado x e y para referirme a las variables y u y v para referirme a los cambios en cada punto y no para lo mismo que en el texto. Es aposta para confundir.

Perdonadme la falta de eficiencia y el probable uso de estructuras ridículas, soy novata.

"""
vectorfield(R, [x1 x2 y1 y2])

Calcula el plano de fase de una ecuación diferencial

Input
- R: matriz de coeficientes
- span: límites del campo [x1 x2 y1 y2]

Output
- x : posición
- y : velocidad
- u : vector de desplazamiento en un punto (x,y)
- v : vector de velocidad en un punto (x,y)
"""
function vectorfield(R, span)
# Span
x = range(span[1],span[2],length=20)
y = range(span[3],span[4],length=20)

# Rejilla de puntos
x = repeat(x', length(x)) # solo repite filas
y = repeat(y, 1, length(y))

# Inicialización
u = zeros(size(x));
v = zeros(size(y));

# Derivada en cada punto
for i = 1:size(x,1)*size(x,2)
  der = R*[x[i]; y[i]];
  u[i] = der[1];
  v[i] = der[2];
end

return (x,y,u,v)
end

Trayectoria: Euler progresivo

Antes de que pintemos el gráfico, en plan bonus voy a contar aquí cómo he simulado el comportamiento del sistema en el tiempo con un Euler progresivo, el método más sencillo para resolver una ecuación diferencial numéricamente.

Básicamente consiste en echar una línea desde el primer punto, que conocemos al ser las condiciones iniciales, hasta un poquito más lejos. La línea en cuestión tendrá la pendiente de la derivada en el punto inicial, que conocemos ya que nos la da la ecuación diferencial. Una vez que hemos llegado al segundo punto, tiramos otra línea y así poco a poco vamos avanzando.

Es decir, en el instante de tiempo i+1 x valdrá lo que valía en el instante i más el cambio:

x_{i+1} = x_i + \dot{x_i}h

siendo h el paso que damos en el tiempo.

He escrito una función que me calcula la posición y velocidad del sistema en el tiempo. ¡Cuidado si le quitáis el amortiguamiento que se vuelve inestable! Temas numéricos eh, no se volvería inestable en la realidad.

"""
euler(R, ic, N, h)

Calcula la respuesta a una ecuación diferencial con Euler progresivo

Input
- R  : matrices en Y(t+1) = Y(t) + h*R*Y(t)
- ic : condiciones iniciales
- N  : número de puntos
- h  : paso

Output
- x, v : desplazamiento y velocidad del sistema en el tiempo
"""
function euler(R, ic, N, h)

Y = zeros(2,N);
Y[:,1] = ic;

for t = 1:(N-1)
 Y[:,t+1] = Y[:,t] + h*R*Y[:,t];
end

x = Y[1,:];
v = Y[2,:];

return (x, v)
end

El dibujo final

Ahora que ya sabemos cuánto cambian tanto la posición y la velocidad en los puntos del plano y hemos calculado la trayectoria para un caso concreto, vamos a pintar. Para ello he usado la función quiver que, dándole la posición y los componentes en ambas direcciones, pinta unas simpáticas flechillas.

using PyPlot

# Datos del sistema
m = 1
k = 1
cc = 2*sqrt(k*m)
xi = 0.3
c = xi*cc

# Datos de la simulación
N = 1000
ti = 0
tf = 50
t = range(ti,tf,length=N)
h = t[2]-t[1]

ic = [1 0]'
R  = [0 1; -k/m -c/m]
F  = zeros(2,N)

# Trayectoria
x,v = euler(R, F, ic, N, h)

clf()
plot(x,v)
gcf()

# Plano de fase
span = [-1 1 -1 1]
x,y,u,v = vectorfield(R, span)
q = quiver(x,y,u,v,color="gray")
xlabel("x")
ylabel("v")
gcf()

Así, nos aparece este gráfico:

phase

Aparte de ser muy cuco, este gráfico nos indica que si le damos cualquier tipo de condición inicial a nuestro sistema, tenderá a pararse: todas las flechitas llevan al punto (0,0). O como diría yo si me quiero hacer la interesante: las raíces de este sistema son complejas con parte real negativa.

Y hasta aquí hemos llegado, os dejo una pregunta de deberes: ¿qué pinta creéis que tendría el plano de fase si le quitásemos el amortiguamiento? Contadme vuestras sospechas.

Referencias

Ejemplos de gráficos en Julia


Os dejo con la canción más bonita del mundo.

Lo que he aprendido: de Matlab a Julia

Hoy vengo a contaros una historia que se inició cuando volví de las vacaciones y descubrí que necesitaba reactivar la licencia de Matlab. Esto llevó a un lío en el que el servicio informático me desinstaló el programa y me instaló una versión hipervieja, en mi departamento me contaron que no, que no tenía que haber contactado con el servicio informático que la licencia la llevaba otra persona, tuve que desinstalar el Matlab viejo, reinstalar el nuevo y activar la licencia como decía esa persona. En definitiva, dos días de trabajo perdidos.

Si ya detestaba Matlab, esto me llevó a odiarlo aun más y buscar una alternativa libre y viable lo antes posible. Ya había andado con GNU Octave antes, pero no quería pasarme a algo que fuera siguiendo a Matlab (y a su sintaxis incoherente, sus funciones cambiantes de versión a versión y demás mierdas variadas que me sacn de quicio). Así que me puse a investigar lo que según el Fediverso es la mejor alternativa libre: el novedosísimo lenguaje de programación Julia cuya versión 1.0 es de agosto de este mismo año.

Sobre el lenguaje

Dicen sus creadores que cuando diseñaron Julia querían un lenguaje abierto (lleva licencia MIT porque se creó en el MIT), tan rápido como C, tan dinámico como Ruby, homoicónico como Lisp, con una notación matemática tan sencilla como la de Matlab, tan generalista como Python, tan útil para cálculo estadístico como R, con tantas capacidades para trabajar con texto como Perl y tan bueno juntando programas entre sí como un intérprete de comandos. Algo sencillo pero que satisficiera a los expertos, que fuera interpretado y compilado. Poquita cosa, vamos.

Lo presentaron en un paper, cosa que a mí ya me gana, y con unos números bastante espectaculares en el benchmark en el que pulen a Matlab en tiempos de computación y se acercan peligrosamente a C++, que toman como referencia. Del lenguaje poco más os puedo decir porque no entiendo nada más, solo que está escrito en una mezcla de C, C++ y Scheme, cosa que lo hace cuanto menos exótico.

En fin, centrándonos en cómo ponernos ya mismo a programar en Julia, no tenemos más que ir a su página de descargas y elegir nuestro método favorito. Yo, de momento, me he instalado el cacharro en el Windows del trabajo.

Diferencias con Matlab

En mi caso concreto que vengo de Matlab, lo que me interesa saber es cuánto de mi conocimiento actual puedo reciclar y la buena nueva es que ¡mucho! La sintaxis de ambos lenguajes es practicamente igual con algunas pequeñas diferencias a las que me he acostumbrado rápido, por ejemplo:

  • Los paréntesis solo se usan para pasar argumentos a las funciones, para elegir elementos de una matriz se usan los corchetes
  • No se copian las matrices, cuando hacemos A = B si cambiamos B cambia A
  • Para aplicar cierta operación elemento a elemento tenemos el operador splat ... o la construcción f.(x)
  • No se agrandan las matrices automáticamente, lo que nos obliga a prealquilar memoria
  • Se pueden tener todas las funciones que se quiera en un mismo archivo

Para los primeros días de locura tenemos hasta un traductor Matlab – Julia que nos ayuda a detectar estas pequeñas diferencias que son difíciles de ver.

Lo que más he notado es que Julia nos obliga a tener buenas prácticas a la hora de programar. Os pongo un ejemplo que me encanta: cuando calculas los autovalores con un método iterativo como Lanczos, te da cada vez cosas diferentes porque empieza con un vector aleatorio si tú, so vago matlabil, no le dices de dónde debe empezar, cosa que deberías hacer porque es tu matriz.

Otro tema chulo y diferente a Matlab es que Julia tiene un gestor de paquetes y debemos decirle explícitamente qué paquetes estamos usando en nuestro script o función con using PAQUETE. De momento he usado Polynomials, Plots, un paquete para cálculo simbólico y otro para matrices dispersas, todos ellos han funcionado muy bien y solo les puedo poner como pega que a veces la documentación es liosa, supongo que porque acaba de salir la versión 1.0. Se pueden además usar paquetes de Python, R o Matlab lo que hace que siempre tengamos recursos para nuestro caso de uso.

Mi entorno juliano actual

Ahora mismo estoy usando la consolilla de Julia (el Julia REPL) en Atom con la extensión Juno. Fui muy reacia a instalar Atom porque consume muchísimos recursos, pero ahora mismo en Windows es la mejor opción. Se instala en nada, tienes todo centralizado y funciona correctamente.

Y tiene un nivel de cuquismo elevado como podéis ver en la siguiente captura de pantalla:

Captura de pantalla de mi entorno para usar Julia en Atom

¿Por qué no desde Emacs?

Creedme cuando digo que lo intenté pero fue imposible. La explicación es sencilla: Windows. La mitad de los paquetes que permiten usar Emacs como IDE de Julia no carrulan en Windows. En teoría con instalar julia-mode y julia-repl o julia-shell-mode debería ser suficiente, pero después de horas infructuosas me rendí. Aun así, mereció la pena la investigación porque descubrí ESS, un modo de Emacs para hacer análisis estadísticos al que espero volver algún día.

Conclusión

Como conclusión os digo que aunque lo he tocado poquito, Julia me encanta, hace lo que yo quiero, es natural programar en él y hay todo tipo de paquetes. Intentaré ir migrando todo mi código de Matlab a Julia e iré compartiendo con vosotros truquillos que vaya descubriendo.

¡En algún momento habrá código que Julia para calcular con elementos finitos gracias a mí! Hasta ese día podéis ir leyendo los enlaces que os dejo en las referencias 🙂

Referencias

Julia: A Fast Dynamic Language for Technical Computing

Julia for Matlab users

Crítica a Matlab en Neuroplausible

How I self-radicalised (an open science manifesto, part 1)

Documentación de Julia

A Deep Introduction to Julia for Data Science and Scientific Computing

Comparación entre Python, Matlab y Julia

Diferencias entre Matlab y Julia

Podcast con los creadores del lenguaje

Chuletilla

Think Julia

Julia Express (pdf)

Lo que he aprendido: leer matrices MMF en Julia

Siguiendo con mi aprendizaje de Julia como sustituto a Matlab, he migrado una función que tenía para leer las matrices de masa y rigidez que exporto de Ansys para luego hacer movidas locas. Otro día os cuento el por qué de la transición (ejem, licencia, ejem) y cómo podéis uniros a mí en el lado oscuro.

Ale, al tema. Primero os voy a contar cómo exportar de Ansys en un par de párrafos, os doy un poco de contexto sobre escribir matrices en archivos de texto y luego vamos al grano.

Exportar matrices de Ansys

Para obtener las matrices de masa y rigidez de Ansys en primer lugar
exportamos el archivo input del modelo con: Tools &gt; Write input file. Luego, a ese input le añadimos antes de la línea /wb,file,end este trocico:

! Gets Stiffness Matrix
! D for double precision real number, C for complex
*SMAT, MatK, D, import, full, file.full, stiff
*export, MatK, mmf, matkMMF.txt ! Exports Stiffness as MMF format

! Gets Mass Matrix
*SMAT, MatM, D, import, full, file.full, mass
*export, MatM, mmf, matmMMF.txt ! Exports Stiffness as MMF format

Seguir leyendo →