Archivo de la categoría: Lo que he aprendido

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: lista de tareas en un documento en LaTeX

Generalmente, tanto cuando escribo algún texto como cuando “programo” me escribo comentarios a mí misma para recordarme lo que tengo que hacer en el futuro. Esos comentarios empiezan por TODO para que me resulte sencillo recuperarlos todos con un simple grep.

Este sistema es simple y eficaz, pero hay ocasiones en las que me gustaría exportar esos comentarios para mí misma al documento que estoy escribiendo, bien para leerlo en su formato final, bien para indicar algunas cosas a mejorar cuando comparto el borrador con algún colaborador.

Como yo siempre escribo en LaTeX, para estos casos he comenzado a utilizar el paquete todonotes que permite anotar el texto en diferentes colores, añadir una señal que indica que falta una figura, o crear una lista de tareas pendientes al final del documento.

Es un paquetillo muy sencillo de utilizar:

  • Se añade una nota a un párrafo en el margen con \todo{Texto de la nota}
  • Se añade una nota en línea con el texto con \todo[inline]{Texto de la nota}
  • Se añade una figura de mentira para avisar de una que falta con \missingfigure{Descripción de la figura}
  • Se incluye una lista de tareas con \listoftodos, esta lista tendrá unos cuadraditos de colores si le pasamos la opción colorinlistoftodos al paquete
  • Soporta varios idiomas, que se mandan como opción cuando llamamos al paquete
  • Si le pasamos la opción disable al paquete nos ventilamos todas las notas

Me parece especialmente práctico si asignamos significado a algunos colores, ya que así podemos ver rápidamente qué nos falta por hacer.

Os dejo para terminar un ejemplo en el que he usado hyperref para al pinchar en la tarea correspondiente nos envíe al lugar donde debemos trabajar:

\documentclass[a4paper]{article}

\usepackage[colorinlistoftodos, spanish, textsize=small]{todonotes}
\usepackage[hidelinks]{hyperref} % hidelinks para que no pinte el cuadrado
\usepackage{lipsum} % texto de prueba

\begin{document}

\lipsum[1]\todo{Reescribir sección}
\lipsum[2]\todo[inline]{Añadir referencias}
\lipsum[3]\todo[inline, color=green]{Comentar los resultados}

\missingfigure{Gráfico mostrando la relación entre las variables}

\listoftodos

\end{document}

Obtendríamos algo así:

tareas

Como siempre, esto es LaTeX, así que podemos cambiar todo lo que nos apetezca. Espero que os resulte útil y que me contéis vuestros trucos.

Referencias

El paquete todonotes en CTAN

How to add todo notes? en StackExchange


Suena:

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: tests en Matlab

He contado más veces que programo a diario sin saber programar. Esto me lleva a situaciones absurdas y unas pérdidas de tiempo sin igual a las que me gustaría poner remedio. Por ese motivo, estoy intentando aprender tácticas que me aligeren el trabajo. Los tests son una de ellas, así que he estado jugando un poco con los tests de Matlab.

Os cuento lo que he hecho y provecho para hacer propaganda de mi campo de estudio: las vibraciones.

Definir la función

Imaginemos que queremos comprobar que una determinada función se comporta de una determinada manera. Voy a usar como ejemplo el MAC (Modal Assurance Criterion), un número que nosotros los ingenieros mecánicos usamos para ver cuánto se parecen dos modos. Dicho así suena muy complejo, pero simplemente es una multiplicación escalar entre dos vectores que previamente hemos normalizado, es decir, que hemos dividido entre su magnitud. Visto en dos dimensiones consiste simplemente en calcular la sombra que un vector proyectará sobre el otro cuando situamos ambos en el mismo origen.

En fin, para dos vectores cualquiera (¡que pueden ser hasta complejos!) se puede calcular el MAC mediante la siguiente función:

function MAC = mac(a,b)
% Calcula el MAC entre los modos a y b
%
% MAC = MAC(a,b)
%
% donde a y b son vectores columna de tamaño Nx1

MAC = abs(a'*b)^2/((a'*a)*(b'*b));

end

Definir los tests

Hemos dicho que el MAC indica cuánto se parecen dos vectores y que representa la sombra de uno sobre el otro. Así, se me ocurren dos manera de comprobar si la función hace lo que debe:

  • Cuando los dos vectores son perpendiculares, el MAC tiene que ser cero porque la longitud de la proyección de uno sobre el otro (la sombra) es cero.

  • Cuando los dos vectores son iguales, el MAC tiene que ser uno porque son lo máximo de parecidos que pueden ser y al ser unitarios la máxima proyección solo puede medir uno.

Ahora toca definir esas condiciones en lenguaje matlabiano. Yo he ido por lo sencillo y he usado vectores bidimensionales:

% TESTS
% Vectores perpendiculares --> MAC = 0
aOrt = [1 0]';
bOrt = [0 1]';

% Vectores iguales --> MAC = 1
aSame = [1 1]';
bSame = [1 1]';

%% Test 1: vectores perpendiculares
assert(mac(aOrt, bOrt) == 0)

%% Test 2: vectores iguales
assert(mac(aSame, bSame) == 1)

El primer trocico son las definiciones, luego viene cada test en una celda de Matlab, un trozo que comienza por dos porcentajes. Este cachito de código lo he guardado como testMac.m y será el archivo que Matlab usará para poner a prueba mi función.

Testar la función

Para lanzar los tests ejecutamos runtest dándole como argumento el ficherillo donde hemos definido los tests:

result = runtests('testMac');

que nos dirá si han fallado los test o no y nos guardará la información que luego podemos ver en una simpática tabla haciendo:

rt = table(result)

Y, hala, ya sé testar. Ahora me falta aprender a definir tests con un mínimo de lógica. Lo próximo será investigar esto mismo en Julia. Os contaré.

Referencias

Testing Frameworks en la documentación de Matlab

Write Script-Based Unit Tests en la documentación de Matlab

Modal space articles de Pete Avitabile

Lo que he aprendido: el buffer de paquetes de Emacs

Años usando Emacs y no ha sido hasta hace unos meses que he descubierto que el buffer de los paquetes, ese al que entramos con M-x list-packages es mucho más fácil de usar de lo que yo pensaba.

Resulta que si te posicionas en la línea de un paquete sin instalar y le das a la i lo marcas para instalar, si vas a uno instalado y le das a la d lo marca para eliminar (delete) y si le das a la U así en general, marca los paquetes que hay que actualizar con una I (upgrade). ¿Y ahora cómo hacemos efectivas estas órdenes? Pues dándole a la x. Así todo lo que queramos borrar se borrará, lo que queramos instalar se instalará y lo que queramos actualizar se actualizará.

En la imagen se ve la I mágica de la que hablo en la columna de la izquierda del todo.

paquetes

En fin, aquí lo dejo, por si hay otro genio como yo por ahí al que le resulta útil. A cuidarse.

Referencias

48.1 The Package Menu Buffer en el manual de Emacs


Solo del punk no se vive, hoy suenan trompetas.

Lo que he aprendido: AUCTeX

Me he hartado de usar más de un programa y he decidido que las pocas veces que escriba en LaTeX a pelo lo haré también en Emacs. Así soy, cada día más simple. Para ello estoy usando AUCTeX, un modo hiperpotente que tiene una manual de solo 130 páginas y que se autodefine como un sofisticado entorno TeX para Emacs.

Lo curioso del tema es que instalé AUCTeX hace más de un año porque fue la única manera de hacer funcionar CDLaTeX, el modo menor que ayuda a crear entornos y movidas matemáticas y que uso con Org. En todo este año no he configurado AUCTeX y he seguido usando Texmaker y Kile por pura inercia. Ahora que lo he puesto a mi gusto y que llevo usándolo una semana, he desinstalado el resto de editores de LaTeX y no creo que vuelva atrás, ¡me encanta!

Unas órdenes básicas

Como siempre, para usar AUCTeX hay que instalarlo y activarlo, a mí se me activa solo al abrir un archivo con extensión .tex, pero si no fuera el caso, M-x LaTeX-mode y adelante. Cuidado con confundirlo con latex-mode, el modo para LaTeX que viene con Emacs. Sí, lo hacemos aposta para liar.

Usar AUCTeX es bastante sencillo, hay un par de atajos de teclado que valen para casi todo y que os listo aquí. La palabreja que pongo en cursiva es con lo que yo relaciono la combinación de teclas para acordarme (que puede o no ser la  idea del autor original).

  • C-c C-c (compile) ejecuta una orden, dando a TAB podemos ver las opciones, hay opciones molonas como Clean all y Clean que eliminan respectivamente todo lo creado por la compilación y solo los archivos auxiliares.
  • C-c C-v (view) muestra el resultado.
  • C-c C-e (environment) introduce un entorno y te va pidiendo los datos que necesite.
  • Si tenemos el activado el parser automático, con C-c C-a (all) deduce qué tiene que usar y compila las veces que haga falta. No viene activado por defecto porque tarda un poco. Esta historieta nos fabrica una carpeta `auto` con cosas dentro, no os asustéis como hice yo. En el archivo de configuración que he puesto más abajo podéis ver cómo se activa.

Como hereje del Emacs que soy os voy a decir algo: no hace falta aprenderse todo esto. Utilizad la GUI y los menús desplegables sin pudor alguno (yo lo hago, un botón tiene un leoncito muy cuqui).

Sigue leyendo

Lo que he aprendido: de SVG a TikZ con Inkscape

Amo Inskape. A pesar de estar en la versión 0.92, funciona maravillosamente, tiene montones de funcionalidades y, por si fuera poco, al más puro estilo software libre, se pueden escribir extensiones que hagan lo que nos apetezca.

Como yo no soy tan hábil como para escribir mis propias extensiones, me aprovecho de la sabiduría de la comunidad para hacer cosas como la que os venía a contar hoy: exportar un svg a TikZ.

Un poco de contexto

Inkscape carrula con imágenes en formato SVG (Scalable Vector Graphics), un formato de gráficos vectoriales definidos en un XML. Que son texto plano, vamos, y podemos abrir una imagen SVG con un editor de texto y ver su definición. La ventaja de usar gráficos vectoriales es que no perdemos calidad al ampliar las imágenes porque los elementos se definen en función de su posición en el dibujo, si el dibujo es más grande las posiciones cambian y los elementos se dibujan de acuerdo a ellas.

TikZ, por su parte, es un paquete de TeX para crear gráficos vectoriales programáticamente1 (qué palabra más maravillosa). De esta manera, los gráficos nos quedan perfectamente integrados en nuestro documento y podemos cambiar el tipo de fuente, tamaño o lo que sea de ambas cosas simultáneamente.

Bien, la cuestión aquí es que usar TikZ es tremendamente difícil y usar Inkscape tremendamente sencillo y ¡LaTeX no tiene soporte para SVG! En general, suelo exportar las imágenes de Inkscape como EPS, manteniendo el texto como LaTeX, pero hay veces que me resulta más práctico tener la definición de la imagen en el propio archivo con el texto.

La extensión svg2tikz

Para esos casos en los que por lo que sea queremos que nuestra imagen esté definida en TikZ, tenemos la extensión de Inkscape svg2tikz. Son tres archivos, dos de extensión de Inkscape y uno de Python que están en la carpeta extensions del repositorio y que tenemos que copiar en la carpeta de extensiones de Inkscape, tal y como indican en las instrucciones. Será o bien /usr/share/inkscape/extensions o ~/.config/inkscape/extensions/2.

Cuidado aquí que yo guardé los archivos directamente del GitHub con Guardar destino como… y luego no me reconocía la extensión. Descargando todo el repo y copiando los archivos correspondientes en cambio sí que me funcionó.

Ahora en el menú Extensiones > Exportar tendremos una nueva opción Export to TikZ path que nos permitirá guardar la imagen como TikZ. Vigilad la ruta porque por defecto no la guarda en la carpeta donde estaba la imagen original sino en la carpeta de las extensiones.

Y ale, ya podemos chulear diciendo que nosotros hacemos todo directamente en LaTeX como hacían los profesores de mecánica de mi universidad que calculaban la base y la ruleta con la fórmula y luego decían que se veían directamente en el dibujo.

Referencias

Chapter 22. Extensions en el manual de Inkscape

svg2tikz

Galería de extensiones de Inkscape

inkscapeMadeEasy


Suena:


  1. Lo usé de tapadillo para pintar estructuras y que me quedaran cuquis lo apuntes, no sé si recordáis. 
  2. Como uno de mis propósitos del 2019 es ir abandonando las tecnologías privativas, voy a intentar hablar lo menos posible de Windows. 

Lo que he aprendido: hourglassing

No sé si os lo he contado alguna vez, pero mi investigación actual implica escribir código de elementos finitos. No, no son elementos muy delgaditos, sino finitos en contraposición a infinitos, es decir que su tamaño es limitado.

Otro día si eso me detengo un poco más en la explicación de los cacharros, que son movidas locas matemáticas, hoy lo que quería es hablar de un problema que se da a menudo en los cálculos de elementos finitos: el hourglassing. Es importante saber a qué se debe y cómo solucionarlo si queremos que nuestra simulación se parezca a la realidad.

Hourglassing: modos sin energía

El hourglassing (relojdearenización) consiste en que aparezcan modos de energía cero al utilizar integración reducida en un elemento. O, dicho en otras palabras, que el elemento sea más blando de la cuenta y se deforme sin que nadie lo toque. De esta manera, cuando calculemos los modos de nuestra estructura, nos aparecerán un montón de modos con un frecuencia cercana a cero. Dependiendo del orden del elemento, habrá más o menos modos de este tipo porque podrá deformarse de diferentes maneras.

Hay que tener cuidado aquí y no pensar que siempre que tenemos muchos modos de frecuencia cero hay hourglassing, ya que puede ser que hayamos metido la pata definiendo las condiciones de contorno, tengamos alguna pieza suelta y sean modos de sólido rígido perfectamente correctos.

Integrando en el elemento

Cuando digo integración reducida, me refiero a usar menos puntos de integración de Gauss de los que nos darían el resultado exacto. Por ejemplo, para un hexaedro lineal, necesitamos dos puntos de integración en cada dirección para obtener una matriz de rigidez exacta, a esto lo llamamos integración completa. Si decidimos usar solo un punto de integración en cada dirección, nuestro cálculo será menos exacto pero más rápido y corremos el riesgo de tener hourglassing. Ocurre lo mismo en el caso cuadrático, con tres puntos en cada dirección tendríamos integración completa y con dos integración reducida.

Podría pensarse que si la integración reducida da esos problemas es absurdo usarla, pero no es así, ya que nos evita un problema que padecen los elementos lineales: el shear locking. Esta otra palabra raruna hace referencia a la aparición de una fuerza cortante de mentira al intentar deformar a flexión un cuadrado o cubo lineal. Como es lineal, sus aristas no pueden curvarse, se comporta como un elemento mucho más rígido y nos da resultados absurdos. Aquí el bailoteo que permite la integración reducida juega a nuestro favor.

Visualizando el hourglassing

Aprovechando mi propio programita, os voy a enseñar la pinta que tienen los resultados de una simulación cuando ha habido hourglassing. En este ejemplo concreto es fácil de ver, ya que la superficie queda como arrugada. También se entiende mejor de dónde viene el nombre.

hourglass

Eso sí, cuidado, ya que no siempre es tan evidente. Lo infalible es calcular los modos: aparecerán los modos reales de la estructura y los espurios de energía cero que no tienen contrapartida física.

Conclusiones

Visto lo visto, os dejo con unos consejos para cuando simuléis, igualmente válidos para los que estéis escribiendo código desde cero y los que uséis software comercial:

  • Si salen modos con valores raros, después de verificar que las condiciones de contorno están bien, yo miraría qué tipo de integración estamos usando.

  • Si cambiando a integración completa los valores siguen sin tener sentido, yo comprobaría el orden del elemento. Si estamos usando elementos lineales con integración completa en un caso de flexión, son demasiado rígidos y van a dar valores de frecuencias naturales muy por encima de los correctos o deformaciones mucho más pequeñas (¡shear locking!). Una opción es pasarse a elementos cuadráticos, aumentando considerablemente el tiempo de cálculo, o considerar los elementos de modos incompatibles, que evitan el shear locking añadiendo restricciones adicionales.

  • Siempre verificaría que los problemas no se deben al tamaño de la malla, no cuesta (casi) nada reducir el mallado y nos quita muchos quebraderos de cabeza. Los elementos chiquiticos son más rígidos y menos propensos al hourglassing.

  • Si todo lo anterior falla, si estuviera en un programa comercial hurgaría en los controles del hourglass o metería amortiguamiento (si estoy en un caso que aplique). Solo como último recurso, todas las mierdas artificiales que le metamos a nuestra simulación tienen tendencia a explotarnos en la cara.

Espero que os resulte útil y poder algún día compartir el código para que juguéis.

Referencias

Why worry about hourglassing in explicit dynamics. Part 1

Why worry about hourglassing in explicit dynamics. Part 2

FEM: Rank deficiency and hourglassing

What is shear locking in FEA?

Lo que he aprendido: ejemplos con código en LaTeX

Estaba currelando en mis traspas para el curso de LaTeX para los doctorandos de mi uni y me pareció buena idea mostrar junto el LaTeX y el resultado para las tablas, ecuaciones y demás mandanga. Sabía que el paquete tcolorbox permite hacer esas historias, pero me parecía demasiado para la poca cosa que quería hacer. Así que hurgué un poco y me encontré con showexpl, un paquetillo que sirve precisamente para hacer lo que yo quería.

No os voy a contar gran cosa porque no deja de ser un añadido a listings que nos da el entorno LTXexample para escribir el código en un lado y el resultado en el otro. De hecho, en las opciones del entorno se puede meter la configuración de los listados y los pies funcionan exactamente igual que los de lstlisting.

Os dejo un ejemplo con los colores de la plantilla de beamer que estoy usando para que juguéis, es bastante claro lo que hace todo:

\documentclass[a4paper]{article}

\usepackage{xcolor}
\usepackage{amsmath}

\usepackage[final]{showexpl}

% Definición de colores
\definecolor{colororange}{HTML}{E65100} % orange
\definecolor{colordgray}{HTML}{795548}  % dark gray for note
\definecolor{colorhgray}{HTML}{212121}  % heavy dark gray for normal text
\definecolor{colorgreen}{HTML}{009688}  % green
\definecolor{colorlgray}{HTML}{FAFAFA}  % background light gray
\definecolor{colorblue}{HTML}{0277BB}   % blue

% Estilo de código
\lstset{
        tabsize=2, 
        backgroundcolor=\color{colorlgray},
        captionpos=b, 
        basicstyle=\small\ttfamily,
        columns=fixed, 
        extendedchars=true, 
        breaklines=true,
        prebreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\hookleftarrow}},
        showtabs=false, 
        showspaces=false,
        keywordstyle=\bfseries\color{colororange},
        commentstyle=\color{colorgreen}, 
        stringstyle=\color{colorblue},
        language={[LaTeX]TeX},
        texcsstyle=*\color{colororange}
}

\begin{document}

% Ejemplo de código
\begin{LTXexample}[pos=r,wide,width=.7,rframe={}]
% - pos: lugar donde saldrá el resultado
% - width: anchura relativa del resultado
% - rframe={}: para que no me ponga un cuadrado alrededor
\begin{equation}
  A = \pi\times R^2
\end{equation}
\end{LTXexample}

\end{document}

que nos crea esto:

A la izquierda el código de LaTeX necesario para escribir una ecuación. A la derecha la ecuación.

Ale, cambiadme parámetros y cosas y me decís qué os sale.

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)