Archivo del Autor: Ondiz

Acerca de Ondiz

Doctora en ingeniería. Loca del software libre. Haciendo GNU/Linux más accesible para mí y de paso para los demás. También hago pan.

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:

En qué ando: mayo

Si abril fue el mes de cuidar la tierra, mayo ha sido el de recoger sus frutos. He comido y regalado lechugas, echado rúcula del huerto a la ensalada y fabricado tortilla de espinaca gigante recién recogida. He plantado tomates, les he puesto palos y me he clavado un pincho del naranjo en la cabeza en el proceso. Todo por ser una agricultora miope. He puesto las calabazas en el jardín y pronto me tocará hacer semillero del resto de las verduras de invierno.

Es una maravilla tener un huerto, no sé cómo no he empezado antes. Es junto con hacer mi propio pan (¡seis años ya de masa madre!) y mis propios yogures (¡desde el 2016!) la mejor decisión que he tomado en la vida. Me falta solo cultivar cereal y ya (casi) podré desconectarme de la sociedad.

La pega es que todo este trajín hortícola no me ha dejado mucho tiempo para dedicarle al blog, solo he escrito sobre ecuaciones diferenciales y notas en LaTeX, y he compartido la receta del bizcocho de calabaza que suelo hacer. Bueno, y me puse intensita un día con la ingeniería. También ando haciendo un estudio sobre los rastreadores en las páginas de los periódicos, pero sabe dios cuándo lo acabaré.

Otra de mis limitaciones blogueriles ha sido que he estado haciendo cosas IRL: me apunté a un curso de comunicación en el currelo; fui a un par de presentaciones de libros; me invitaron a un congreso; estuve hablando con la gente de Finantzaz haratago y les di la turra para que se pasasen a las redes libres (y se pasaron), descubrí por mi hermano un proyecto de un internet paralelo basado en pollos… Vamos, lo normal.

Algunas cosas interesantes

Artículos

Enlazo aquí algunas cosas que he leído, hilando con lo que comentaba antes del rastreo recomiendo usar uBlock Origin y Privacy Badger para entrar en los artículos, sobre todo los de los periódicos.

Libros

Tras mi incursión en las letras he vuelto a las matemáticas, sigo poco a poco con Gödel, Escher, Bach que leo siguiendo el curso del MIT al respecto y estoy leyendo a desgana Weapons of math destruction, que me recomendaron. Como académica elitista que soy, detesto los títulos sensacionalistas y la superficialidad en los análisis.

Sobre todo y ante todo me he obsesionado con un libro sobre la autosuficiencia: The new complete book of self-suffienciency. The classic guide for realists and dreamers y me dedico a leer sobre verduras, cerveza y fabricar mi propio horno con mis propios ladrillos. He llegado a la conclusión de que es tarde para volverme normal así que voy a abrazar mi locura.

Venga, va, lo confieso: no he abandonado las letras. Siguiendo con mi vieja costumbre de leer cosas que no entiendo, he leído The rhetorics of economics y ando enamorada de Contar es escuchar, una recopilación de textos de no ficción de Ursula K. Le Guin en el que se reflexiona sobre la lectura y la escritura. Son especialmente hermosas las partes que tratan la relación entre el escritor y lector, la escritora y la lectora en este caso. También cometí el grave error de entrar en el Diccionario panhispánico de dudas; ahí sigo.

Podcasts


Ya es casi verano, bailemos.

De mayor quiero ser ingeniera

UPDATE: thanks to JordiGH, there is a translation of the post to English after the original


Cuando era pequeña leía y dibujaba todo el rato. Leía tanto que me tuvieron que pasar a la sección de adultos antes de la cuenta porque me había quedado sin opciones. Y dibujaba. El manga se puso de moda cuando yo tenía unos once o doce años así que aprendí a dibujar por mi cuenta y me entretenía en clase. Prestar atención no merecía la pena y hablar con mis compañeros… es largo de contar.

A pesar de esto (o precisamente por esto), siempre me han fascinado la ciencia y, sobre todo, las matemáticas. Hay algo hermoso en poder plasmar la realidad en unos símbolos en un papel, en predecir el comportamiento de un objeto, en comprender cómo funciona una máquina. Así, como era lista y las letras no valen para nada (¿a quién no le han dicho esto?) acabé estudiando una ingeniería.

De la carrera mejor no hablar: me dieron soluciones a problemas que no tenía, los problemas que tenía se quedaron sin resolver. Pero la saqué. Hice un master. Una tesis. Un postdoc. Los problemas siguieron ahí. La rueda siguió adelante.

Durante años pensé que no era una ingeniera de verdad porque seguía dedicándole mi tiempo libre al arte, a la literatura, no como todas esas personas a mi alrededor que hablaban de coches y estructuras y código en su tiempo libre, que leían sobre nuevos avances tecnológicos para no quedarse atrás mientras yo lo único que quería era salir a la calle con mis acuarelas. Pensé que no podía ser una ingeniera de verdad porque pensaba que el progreso por el progreso no nos llevaba a ningún sitio, porque quería desarrollar tecnología humana y respetuosa con el planeta, porque no quería vivir en un mundo controlado por un puñado de empresas, porque quería resolver problemas reales que afectasen a personas reales. Porque quería, en definitiva, hacer del mundo un lugar mejor. 

Era una loca por no querer un trabajo bien pagado, a quién le importaba que fuera repetitivo e inútil.

En las letras no me trataron mejor: era el enemigo. Destruía empleos. Era cabezacuadrada y carente de creatividad. Quería someter al ser humano a los designios de la tecnología. Nunca quise, pero no importó. Y me quedé en el gris: demasiado artística para encajar en la ingeniería, demasiado ingenieril para encajar en las artes o en las letras.

De esta situación me sacó, como tantas veces en la vida, la lectura. La ingeniería como disciplina donde señores de traje discuten sobre hojas de cálculo es un invento reciente. Hace años había hombres y mujeres que creaban cosas a partir de su imaginación, que escribían textos hermosos y compartían sus conocimientos. Personas valoraban la estética y la elegancia de sus inventos y no eran por ello peores científicas o tecnólogas. Personas que cambiaron de verdad el mundo.

Así, tantos años después de elegir una carrera prácticamente al azar, después de haber pensado que me había equivocado de vida, ahora digo, por fin convencida: de mayor quiero ser ingeniera.


When I was little, I used to read and draw all the time. I read so much that they had to move me to the adult section prematurely because I had ran out of options. And I drew. Manga got popular when I was around eleven or twelve, so I learned to draw on my own and I amused myself in class. It wasn’t worth it to pay attention, and talking with my classmates… it’s a long story.

Despite this (or precisely due to this), I’ve always been fascinated by science, and especially by mathematics. There is something beautiful about being able to capture reality with some symbols on paper, about predicting the behaviour of an object, about learning how a machine works. Thus, since I was smart and art is worthless (who hasn’t heard this?), I ended up studying engineering.

Best to not talk about my undergrad: they gave me solutions to problems I didn’t have, and the problems I did have remained unsolved. But I did it. I did a master’s. A thesis. A postdoc. The problems were still there. The wheel moved forward.

For years I thought that I wasn’t a real engineer because I kept spending my free time on art, on literature, not like all the people around me who talked about cars and structures and code in their free time, who read about new technological advances so they wouldn’t be left behind while all I wanted was to go out to the street with my watercolours. I thought I couldn’t be a real engineer because I thought that progress for its own sake would take us nowhere, because I wanted to develop a humane technology that would respect the planet, because I didn’t want to live in a world controlled by a handful of corporations, because I wanted to solve real problems that impacted real people. Because I wanted, unequivocally, to make the world a better place.

I was crazy for not wanting a well-paying job. Who cares if it was repetitive or useless.

In the humanities they were not kinder to me: I was the enemy. I destroyed jobs. I was a square who lacked creativity. I wanted to make human beings submit to the whims of technology. I never wanted to, but it didn’t matter. And I remained in a grey area: too artistic to fit into engineering, too engineery to fit into art or literature.

As so often it happens in life, reading brought me out of this situation. Engineering as a discipline in which men in suits discuss spreadsheets is a recent invention. Some years ago there were men and women who created things out of sheer imagination, who wrote beautiful texts and shared their knowledge. People valued aesthetics and the elegance of their inventions, and this didn’t make them worse scientists or technologists. People who really changed the world.

And so, after picking a major practically at random, after thinking that I had chosen the wrong life, I say, at last with conviction: I want to be an engineer when I grow up.

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.

Camino a la panificación: bizcocho de calabaza

Hoy traigo un bizcocho que suelo hacer a menudo, es poco dulce, húmedo y denso, como me gustan a mí los bizcochos. Suelo aprovechar para asar una calabacilla de campana cuando hago pan, se aprovecha mejor el calor y luego hay deliciosa calabaza asada para bizcochar, risottear o hacer delicioso hummus calabácido (con la tortilla de patata y el miso mi alimento favorito del mundo).

El último que hice tuvo el mérito añadido de proceder de una calabaza producida por mi huerto. Metros y metros cuadrados de plantas de calabaza que produjeron exactamente dos calabazas. En fin, os dejo mi proceso de crear bizcocho naranja.

Bizcocho de calabaza

  • Dificultad: echar productos en un bol
  • Imprimir

Verdura en el dulce




Ingredientes

  • 280g calabaza asada
  • 200g yogur líquido casero
  • 60g aceite
  • 2 huevos
  • 260g harina integral de espelta
  • 100g azúcar
  • 1 cucharadita levadura
  • 1 cucharadita bicarbonato
  • 2 cucharaditas especias (canela, jengibre, nuez moscada)

Direcciones

  1. Bato los huevos con el azúcar con brío con las varillas.
  2. Añado el aceite y el yogur y bato más.
  3. Agrego la harina con la levadura, el bicarbonato y las especias y mezclo, esta vez con una cuchara, lengua o similar.
  4. Uno la pulpa de calabaza asada a la mezcla con la cuchara.
  5. Horneo en horno precalentado a 180C durante 40′.

¡La próxima vez que lo haga a ver si le saco foto!

En qué ando: abril

El mes de abril ha estado marcado por mi trabajo en el huerto: plantar, desherbar, regar y ver crecer la vida. Abril ha sido un mes de introspección, así que he escrito poco y leído mucho. Solté un poco de lastre escribiendo una cosa rara sobre la vida gracias a la cual me ha escrito mucha gente. Es curioso cómo ha llenado la gente los huecos y ha interpretado el texto según su propia experiencia. Es lo que quería y me alegro de haberlo conseguido.

Lechuga de hoja de roble en primer plano, rúcula e hinojo detrás

Luego ya volví al Emacs, con un truquillo sobre el buffer de paquetes y aprendí a hacer tests unitarios en Matlab.

Algunas cosas interesantes

Este mes he leído muchísimo. He recuperado el vicio lector de cuando era adolescente después de unos cuantos años de desconexión. Supongo que tiene relación con leer cosas buenas. Listo aquí unas pocas.

Sigue leyendo

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.

El mundo gira gira

Uno de los motivos por los que empecé este blog es porque tengo una necesidad física de escribir. Esta entrada obedece a esa necesidad.


1

Siento a menudo que el mundo va demasiado rápido para mí. Me ocurrió por primera vez con catorce años al volver de las vacaciones de verano. Yo seguía siendo yo mientras que la vida en el colegio había avanzado a un lugar temporal inalcanzable para mí, un lugar repleto de códigos que yo no entendía y normas absurdas que me negaba a respetar, que todavía me niego a respetar. Conseguí más tarde seguir más o menos el ritmo que me imponía la realidad, nunca del todo, manteniendo siempre un pie fuera del torbellino, respetándome, respetando mis tiempos.

Ahora el mundo ha llegado a un nivel de aceleración que, como hace media vida, otra vez no soy capaz de comprender. Veo a la gente correr a mi alrededor para mantenerse en el mismo lugar. Trabajar, trabajar, trabajar. Sin pararse un momento a pensar, preguntarse hacia dónde corren o qué es eso que buscan con tanta ansia.

Por esto reivindico la lentitud, disfrutar del camino. Aunque nunca llegue al destino. ¿Hay acaso un destino?

2

La consecuencia de ser demasiado sensible es preocuparse de todo todo el rato. Me imagino accidentes de todo tipo cuando conduzco, pienso que me van a mandar callar, decirme que yo qué sé, que no tengo derecho a estar ahí, que quién me he creído yo que soy. Me imagino todos los escenarios posibles cuando voy a hacer cualquier cosa con diferentes variaciones de cada uno de ellos en un árbol gigante de universos paralelos que responden a la pregunta ¿y si? hasta el infinito. Pienso que voy a morirme porque me late el corazón a mil por hora en medio de la noche pero es solo un ataque de pánico. Y se acaba. Y no pasa nada.

Mi truco después de tantos años es ponerme en el peor de los casos y ver que no es para tanto. Que sobreviviré. Que he sobrevivido.

3

Lo que me mantiene cuerda es tener la cabeza en las nubes pero los pies en la tierra. Siempre digo que el arte nos va a salvar, estoy aquí ahora por la pintura, la literatura, la música. Pero también por el pan y el huerto, ver crecer la vida me hace estar atada al tiempo y al espacio. Cultivar la mente y cultivar la tierra. Y crear lazos humanos, conectarnos, escucharnos, cuidarnos. No hay nada más revolucionario.

Hay que seguir empujando la piedra a pesar de que vaya a volver a caer por la ladera, precisamente porque va a volver a caer por ladera.

Así que cuando tengáis miedo, cuando penséis que no hay luz, que el pozo es demasiado profundo, venid a sentaros conmigo en el borde del precipicio. Veremos amanecer.

Aquí siempre hay sitio.