Análisis Topológico de Datos y Aprendizaje Automático en la Detección De Ríos Atmosféricos

Análisis Topológico de Datos y Aprendizaje Automático en la Detección De Ríos Atmosféricos David Vázquez Fernández Análisis Topológico de Datos y Aprendizaje Automático en la Detección De Ríos Atmosféricos David Vázquez Fernández Memoria presentada como parte de los requisitos para la obtención del título de Máster Universitario en Matemáticas por la Universidad de Sevilla. Tutorizada por Prof. María del Rocío González Díaz Índice general Resumen 1 Abstract 3 1. Introducción 5 2. Conceptos previos 9 2.1. Homología persistente en complejos cúbicos . . . . . . . . . . . . 9 2.1.1. Complejos cúbicos . . . . . . . . . . . . . . . . . . . . . . 9 2.1.2. Homología en complejos cúbicos . . . . . . . . . . . . . . . 12 2.1.3. Complejos cúbicos filtrados a partir de una matriz de datos (construcción T) ....................... 16 2.1.4. Homología persistente en complejos cúbicos filtrados . . . 19 2.2. Imágenes de persistencia . . . . . . . . . . . . . . . . . . . . . . . 22 2.3. Modelos de aprendizaje automático. Clasificación supervisada . . 26 2.3.1. Sobreajuste, división de conjuntos, validación cruzada. . . 28 2.3.2. Elección del mejor método de aprendizaje automático con validación cruzada anidada LTO . . . . . . . . . . . . . . 29 2.3.3. Evaluación de la calidad de una clasificación supervisada binaria............................. 30 ii análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos 2.4. Máquina Vector Soporte . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.1. Datos separables linealmente . . . . . . . . . . . . . . . . . 32 2.4.2. Datos cuasi-separables linealmente . . . . . . . . . . . . . 33 2.4.3. Datos no linealmente separables . . . . . . . . . . . . . . . 37 2.5. KNN, método de los K Vecinos más Cercanos . . . . . . . . . . . 38 3. Estado del arte 41 3.1. Datosempleados ........................... 41 3.2. Esquemadelmétodo ......................... 43 3.3. Primera etapa, análisis topológico de datos . . . . . . . . . . . . . 44 3.3.1. Gráficos de evolución . . . . . . . . . . . . . . . . . . . . . 46 3.3.2. Ejemplo simple del funcionamiento del algoritmo . . . . . 49 3.4. Segunda etapa, clasificación con máquina vector soporte . . . . . 51 3.5. Resultados............................... 54 4. Método propuesto 57 4.1. Datosempleados ........................... 58 4.2. Esquema del método y diferencias con el método existente . . . . 60 4.3. Primera etapa, análisis topológico de datos . . . . . . . . . . . . . 61 4.3.1. Obtención de función de filtrado . . . . . . . . . . . . . . . 63 4.3.2. Ejemplo simple del funcionamiento del algoritmo . . . . . 66 4.3.3. Generación de vectores resultado. Imágenes de persistencia 68 4.3.4. Esquema final y ejemplo de la primera etapa . . . . . . . . 71 4.4. Segunda etapa, método de los K Vecinos más Cercanos . . . . . . 73 4.5. Implementación y código . . . . . . . . . . . . . . . . . . . . . . . 77 ÍNDICE GENERAL iii 4.6. Resultados............................... 77 4.6.1. Importancia de las dimensiones empleadas en homología y del método de aprendizaje automático escogido . . . . . . 78 5. Comparación de métodos: Estabilidad y resultados numéricos 81 5.1. Estabilidad .............................. 81 5.1.1. Estabilidad de la primera etapa del método propuesto . . . 83 5.1.2. Carencia de estabilidad del método de Muszynski, G. et al. 91 5.2. Comparación de resultados numéricos . . . . . . . . . . . . . . . . 95 6. Conclusiones y trabajo futuro 99 6.1. Conclusiones.............................. 99 6.2. Trabajofuturo ............................ 100 Bibliografía 103 Resumen Un río atmosférico (RA) es una banda alargada y estrecha de humedad muy concentrada en la parte inferior de la troposfera que va desde los trópicos hasta un continente. La selección de umbrales para definir cuantitativamente un RA continúa siendo un problema abierto. Es por ello que se han hecho esfuerzos en desarrollar métodos de detección automática de RA en datos climáticos que no dependan de umbrales prefijados. En este trabajo, primero, estudiamos un método propuesto que extrae información topológica sobre componentes conexas en imágenes climáticas y luego, emplea una Máquina Vector Soporte para clasificar esa información según se corresponda o no con que haya un RA en la imagen de la que proviene. Después, proponemos una mejora del método anterior, extrayendo información topológica, no solo sobre componentes conexas, sino sobre huecos, y utilizando el método de los K Vecinos más Cercanos para clasificar esa información. Ambos métodos no fijan previamente umbrales. Posteriormente, demostramos que el método propuesto en este trabajo es estable considerando un paso suficientemente pequeño y que el método ya existente no lo es. Finalmente, al comparar los resultados numéricos de los dos métodos, se concluye que nuestro método es el mejor en la mayoría de medidas de bondad de clasificación calculadas. El código de implementación del método desarrollado en esta memoria está disponible en: https://github.com/davafe/Metodo-deteccion-rios-atmosfericos 8análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Podemos destacar dos métodos propuestos recientemente, que, con otros fines, emplean técnicas parecidas a las del método que proponemos nosotros. En [10], se usa homología persistente (en dimensión 1) en complejos cúbicos construidos a partir de una imagen en escala de grises, para estudiar ciclos diarios en la nubosidad de los ciclones tropicales. En [11], se propone un método para clasificar imágenes en escala de grises que, primero, emplea homología persistente en complejos cúbicos y luego aprendizaje automático. En la última parte de este trabajo, se comparan los dos métodos, el analizado y el propuesto, en términos de estabilidad y resultados numéricos. Se demuestra que el método propuesto en este trabajo es estable considerando un paso suficientemente pequeño y que el método ya existente no lo es. Además, al comparar los resultados numéricos de los dos métodos, se concluye que nuestro método es el mejor en la mayoría de medidas de bondad de clasificación calculadas. Por tanto, el método que proponemos en este trabajo, es mejor en términos de estabilidad y resultados numéricos que el método ya existente. A continuación, detallamos el contenido de cada capítulo de la memoria. En el capítulo 2 se introducen los conceptos previos necesarios para el trabajo y en el capítulo 3 se describe y analiza el método desarrollado en [7]. Posteriormente, en el capítulo 4 se explica el nuevo método propuesto y en el capítulo 5 se comparan la estabilidad y resultados numéricos de ambos métodos. 2 Conceptos previos En este capítulo explicamos conceptos previos que se han empleado en este trabajo. En la sección 2.1 se explican los conceptos previos necesarios para calcular la homología persistente en complejos cúbicos, basándonos en [10], [12], [13] y [14]. En la sección 2.2 se definen las imágenes de persistencia y se comenta su utilidad, basándonos en [15]. En la sección 2.3 se explican conceptos sobre modelos de aprendizaje automático, en particular de clasificación supervisada. Nos hemos basado en [17] y [18]. Finalmente, se exponen los métodos de aprendizaje automático que se usan: Máquina Vector Soporte (en la sección 2.4) y los K Vecinos más Cercanos (en la sección 2.5). Para introducir la Máquina Vector Soporte nos hemos basado en [19] y para introducir el concepto de los K Vecinos más Cercanos en [18]. 2.1 Homología persistente en complejos cúbicos En la subsección 2.1.1, introducimos los complejos cúbicos y en la subsección 2.1.2 explicamos la homología en complejos cúbicos. En la subsección 2.1.3 explicamos cómo podemos construir complejos cúbicos filtrados a partir de una imagen en escala de grises y, en la subsección 2.1, vemos cómo se puede aplicar homología persistente en esos complejos cúbicos filtrados. 2.1.1 Complejos cúbicos Definición 2.1. Un intervalo elemental es un intervalo cerrado I⊂Rde la forma I= [l, l + 1] oI= [l, l]para cierto l∈Z. Para simplificar la notación escribimos [l]=[l, l]en los intervalos que contengan un solo punto. 10 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Los intervalos elementales que consisten en un solo punto son degenerados, mientras que los de longitud 1son no degenerados. Definición 2.2. Un cubo elemental Qes el producto de intervalos elementales, i.e., Q=I1×I2× · · · × Id⊂Rddonde Iies un intervalo elemental para todo i∈ {1, . . . , d}. El conjunto de todos los cubos elementales en Rdse denota por Kd y el conjunto de todos los cubos se denota por K, i.e., K=∪∞ d=1Kd. Definición 2.3. Sea Q=I1×I2× · · · × Id⊂Rdun cubo elemental. El número de embebimiento de Q, denotado por emb Q, es la dimensión del espacio euclidiano en el que está contenido, en este caso emb Q=dya que Q⊂Rd. El intervalo Iies referido como la componente i-ésima de Qy se escribe Ii(Q). La dimensión de Qse define como el número de componentes no degeneradas de Qy se denota por dim Q. Observar que si emb Q=dentonces Q∈ Kd. Además, definimos Kk:= {Q∈ K| dim Q=k} y Kd k:= Kk∩ Kd para todo k, d ∈ {0,1, . . . } Observación 2.1.Claramente se tiene que, para cualquier cubo elemental Q: 0≤dim Q≤emb Q. En la figura 2.1 podemos ver ejemplos de cubos elementales con distintos números de embebimiento y dimensión. Definición 2.4. Sean Q, P ∈ K. Qes una cara de Psi Q⊂P. Esto se denota como Q⪯P. Qes una cara propia de Psi Q⪯PyQ=P, lo que se denota Q≺P. Qes una cara primaria si Q⊂Pydim Q= dim P−1. Ejemplo 2.1.Sea Q= [1,2] ×[1],R= [1] ×[1] yP= [1,2] ×[1,2]. Entonces R⪯Q⪯P⪯P,R≺Q≺PyQes una cara primaria de PyRes una cara primaria de Q.Rno es una cara primaria de P. 2. conceptos previos 11 Figura 2.1: Ejemplos gráficos de cubos en RyR2. El cubo de arriba tiene número de embebimiento y dimensión 1, los cubos del medio tienen número de embebimiento 2pero dimensión 1y, el cubo de abajo tiene dimensión y número de embebimiento 2. Imagen sacada de [12]. 12 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Definición 2.5. Un conjunto X⊂Rdes cúbico si Xpuede ser escrito como una unión finita de cubos elementales. Dado X⊂Rdconjunto cúbico entonces adoptamos la siguiente notación: K(X) := {Q∈ K|Q⊂X} y Kk(X) := {Q∈ K(X)|dim Q=k}. AK(X)se le llama complejo cúbico asociado a X. Un L⊂ K(X)es un subcomplejo cúbico de K(X)si es un complejo cúbico, i.e., si L=K(Y)para un Y⊂X que sea conjunto cúbico. De manera análoga a los grafos o a los complejos simpliciales, a los elementos de K0(X)se les llama vértices de X, a los elementos de K1(X)se les llama aristas de Xy, de manera general, a los elementos de Kk(X)se les llama k-cubos de X. Ejemplo 2.2.Sea X= [0,1] ×[0,1] ∈R2.Xes un cubo elemental y, por tanto, es un conjunto cúbico. Es fácil de comprobar que: K2(X) = {X}, K1(X) = {[0,1] ×[1],[0,1] ×[0],[0] ×[0,1],[1] ×[0,1]}, K0(X) = {[0] ×[0],[0] ×[1],[1] ×[0],[1] ×[1]}, K(X) = K2(X)∪ K1(X)∪ K0(X). Observación 2.2.El análogo a los complejos cúbicos si trabajamos con símplices es el complejo simplicial, donde el conjunto cúbico es lo análogo al espacio subyacente de los complejos simpliciales. Cuando se calcula la homología simplicial, realmente se está calculando la homología en función de la topología del espacio subyacente; de manera análoga, cuando se calcula la homología cúbica, se hace en función de la topología (como subespacio de Rdcon la topología usual) de los conjuntos cúbicos. 2.1.2 Homología en complejos cúbicos La homología es una herramienta estándar en topología algebraica que nos da un grupo para cada dimensión k∈ {0,1, . . . }y para un espacio topológico X 2. conceptos previos 13 dado. Cada una de las dimensiones mide distintas propiedades del espacio. En el caso particular de este trabajo, estaremos interesados solamente en homología 0 y1-dimensional, ya que estamos trabajando con datos bidimensionales (imágenes en escala de grises). Definamos la homología para un conjunto cúbico. La homología depende del grupo Gque se considere para los coeficientes. En este trabajo vamos a considerar G=Z2y así no vamos a tener en cuenta la orientación. Así, las definiciones que vamos a hacer a continuación asumimos que el grupo que se considera para los coeficientes es Z2.1 Definición 2.6. Sea Xun conjunto cúbico, dado k∈ {0,1, . . . }, una k-cadena ccon coeficientes en Z2es una combinación lineal formal de k-cubos en X, i.e., c=X Qi∈Kk(L) aiQi con coeficientes ai∈Z2para todo Qi∈ Kk(L). Podemos sumar estas cadenas de la siguiente forma X Qi∈Kk(X) aiQi+X Qi∈Kk(X) biQi=X Qi∈Kk(X) (ai+bi)Qi. De esta forma, Ck(X) := {PQi∈Kk(X)aiQi|ai∈Z2∀Qi∈ Ki(X)}es el conjunto formado por todas las k-cadenas con coeficientes en Z2formadas por combinación lineal de k-cubos de X. Se puede comprobar que Ck(X)es un grupo abeliano libre con base Kk(X). Dado c∈Ck(X), diremos que la cadena ces de dimensión k. Dado c=PQi∈Kk(X)aiQi∈Ck(X), el soporte de ces |c|:= ∪{Qi∈ Kk(X)|ai= 0}. En la definición anterior estamos denotando de igual manera a los k-cubos que a un objeto algebraico que los identifica. Una vez tenemos una cadena, una combinación lineal de objetos algebraicos asociados a k-cubos, obtenemos su soporte, el espacio topológico formado por 1Si se quiere ver de manera mucho más general la homología cúbica con coeficientes enteros y con resultados y demostraciones, consultar [12]. La adaptación al caso de Z2está basado en la adaptación hecha en [10]. 14 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos la unión de los k-cubos cuyos objetos algebraicos asociados tienen coeficiente no nulo en la cadena. Como los coeficientes de las cadenas pertenecen a Z2, solo pueden valer o 0o 1, esto conlleva a que no hay dos cadenas con el mismo soporte, i.e., a la hora de combinar algebraicamente cubos solo estamos teniendo en cuenta los cubos que combinamos y por lo tanto no tenemos en cuenta orientación, multiplicidad, etc. de los cubos. Ejemplo 2.3.Consideremos X= [0,1] ×[0,1], un conjunto cúbico. Hemos visto en el ejemplo 2.2 que las aristas de Xson K1(X) = {A1:= [0,1] ×[1], A2:= [0,1] ×[0], A3:= [0] ×[0,1], A4:= [1] ×[0,1]}. Un ejemplo de cadena 1dimensional c∈Ci(X)sería: c=A1+A2+A3+A4, dónde |c|es el borde del cuadrado X.X(como objeto algebraico) no deja de ser una cadena 2dimensional. En la siguiente definición definimos una aplicación que a cada cadena le asigna una cadena “borde”, extendiendo los bordes topológicos de los soportes de las cadenas a las propias cadenas. Definición 2.7. Dado k∈ {1,2, . . . }y un conjunto cúbico X, definimos el k-operador borde para el conjunto cúbico X: δk:Ck(X)−→ Ck−1(X) Es un homomorfismo de grupos abelianos libres, que se define para los elementos Q∈ Kk(X)y se extiende por linealidad al resto de Ck(X). Dado Q∈ Kk(X), definimos el borde de Qcomo la suma de sus caras primarias δk(Q) := X P≺Q, dim P=dim Q−1 P . Por tanto, dado ∈Ck(X),δk(c) = PQi∈Kk(X)aiδk(Qi). Esta última definición está particularizada para el caso de considerar coeficientes en Z2. Con esta proposición se ve, de manera algebraica, que el borde de un borde de una cadena es 0. 2. conceptos previos 15 Proposición 2.1.Sea Xun conjunto cúbico y k∈ {1,2, . . . }, entonces δk−1◦δk= 0 . Podemos encontrar una demostración de esta proposición (en el caso más general de coeficientes enteros) en [12]. Ahora ya estamos preparados para definir los grupos de homología. Definición 2.8. Sea Xun conjunto cúbico y k∈ {0,1, . . . }, el grupo de homología de dimensión kde Xes: Hk(X) := Ker(δk)/Im(δk+1). Definimos el grupo de los k-ciclos como Zk(X) := Ker(δk)y el grupo de los kbordes como Bk(X) := Im(δk+1). Observemos que los grupos de homología están bien definidos, ya que, por la proposición 2.1, tenemos que Im(δk+1)⊂Ker(δk). Denotaremos indistintamente a Hk(X), como Hk(X)(grupo de homología de dimensión kdel conjunto cúbico X) y como Hk(K(X)) (grupo de homología de dimensión kdel complejo cúbico K(X)). Esto está en línea con lo comentado en la observación 2.2. Un k-ciclo de Xes un elemento cuyo borde es nulo y un k-borde de Xes un elemento que es borde de una cadena k+ 1 dimensional. Así, las clases no nulas de Hk(X)tienen como representantes k-ciclos que no son k-bordes. Con el siguiente teorema, demostrado en [12], vemos qué información nos aporta H0. Teorema 2.1. Sea Xun conjunto cúbico, entonces H0(X)es un grupo abeliano libre. Además, si {Pi∈ K0(X)|i∈ {1, . . . , n}} es una colección de vértices en Xtal que hay un vértice por cada componente conexa de X, entonces: {[Pi]∈H0(X)|i∈ {1, . . . , n}} es una base de H0(X).[Pi]es la clase de H0(X)en la que está Pi Con este teorema podemos interpretar que existe un elemento de la base de H0(X)por cada componente conexa de X, así que H0(X)nos da información sobre el número y cuales son las componentes conexas de X. 16 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos La interpretación que se le puede dar a H1(X)es similar: el núcleo de δ1, Z1(X), está generado por 1-cadenas cuyo borde sea nulo pero no sean borde de una 2-cadena. Es decir, H1(X)está generado por cadenas que representan uniones de aristas cuyo borde sea nulo (i.e., estas uniones tienen los vértices repetidos un número par de veces, por lo tanto es una unión de aristas que es cerrada) y que no sean borde de 2-cadenas (i.e., que no sean bordes de una unión de 2-cubos, por lo tanto que esa unión de aristas cerrada no esté rellenada por 2-cubos). Por tanto un elemento no nulo de H1(X)es una combinación de huecos en Xy de la misma manera, hay un elemento de la base de H1(X)por cada hueco en X. 2.1.3 Complejos cúbicos filtrados a partir de una matriz de datos (construcción T ) En esta subsección veremos lo que es un complejo cúbico filtrado y como podemos construir un complejo cúbico a partir de una matriz de datos. Una matriz de datos se puede ver como una imagen digital en escala de grises de dimensión 2, dónde el “color” de cada píxel (no deja de ser un número) está determinado por el valor de los elementos de la matriz. Para ello veremos cómo se construyen de manera general complejos cúbicos filtrados a partir de una imagen en escala de píxeles de cualquier dimensión. Con esta construcción podremos aplicar la homología persistente (que explicaremos en la subsección 2.1.4) a una matriz de datos para captar los patrones topológicos que hay en esa matriz de datos. Empezamos viendo lo que es un complejo cúbico filtrado. Definición 2.9. Sea X⊂Rqun conjunto cúbico, una función f:K(X)−→ R es monótona si P⪯Qimplica f(P)≤f(Q). Definición 2.10. Un complejo cúbico filtrado es un par (K(X), f)formado por un complejo cúbico asociado a un conjunto cúbico X⊂Rq,K(X), y una función monótona f:K(X)→R. A fse le llama función de filtrado. Observación 2.3.Observar que la condición de monotonía implica que para todo r∈R,f−1(−∞, r]incluya a todos los cubos de dimensión qcuya imagen por f es menor que ry a todas las caras propias de estos cubos (ya que su imagen tiene que ser menor o igual a la de los cubos por monotonía). Esto implica que 2. conceptos previos 17 f−1(−∞, r]es un subcomplejo cúbico, en particular, f−1(−∞, r] = Kq [ k=0 [ P∈Kk(X) f(P)≤r P. Veamos qué es una imagen dimensional en escala de grises. Definición 2.11. Una imagen (digital) d-dimensional en escala de grises de tamaño (n1, n2, . . . , nd)es un array con valores en R, i.e., I ∈ Mn1n2...nd(R). De manera equivalente, es una función definida en una rejilla: I:I=J1, n1K×J1, n2K× · · · × J1, ndK−→ R donde, dado i∈ {1, . . . , d}, J1, niK:= {k∈N|k∈[1, ni]}. El conjunto de índices Ise llama el dominio de la imagen I. Si d= 2, los elementos de Ise les llama píxeles y, si d≥3, se les llama vóxeles. Es claro que una matriz n×mde datos es una imagen 2-dimensional en escala de grises de tamaño (n, m). Para inducir una topología que represente una imagen en escala de grises y represente la conexión entre sus vóxeles/píxeles no podemos quedarnos con la topología relativa en I⊂Zd⊂Rd, ya que con esa topología la imagen es totalmente disconexa. Es por ello que vamos a inducir la topología en la imagen a través de complejos cúbicos filtrados. Existen dos formas de construir un complejo cúbico filtrado a partir de una imagen ddimensional. Una es la construcción V donde cada vóxel/píxel representa un vértice del complejo cúbico y la función se define en los vértices según el valor del píxel correspondiente y en los cubos de dimensión superior como el máximo de la función en los vértices que son cara de esos cubos. Sin embargo, en este trabajo vamos a usar la construcción T, dónde cada vóxel/píxel representa un d-cubo del complejo cúbico y la función se define en los d-cubos según el valor del píxel correspondiente y en los cubos de dimensión inferior a dcomo el mínimo de la función en los d-cubos que los tienen como cara propia. Podemos ver un ejemplo de comparativo de las dos construcciones en la figura 2.2 24 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Definición 2.17. Sea Bel conjunto formado por todos los puntos de un diagrama de persistencia que tienen multiplicidad positiva y no tienen persistencia infinita, sea g:R2−→ R una función diferenciable a trozos tal que g(x, 0) = 0 para todo x∈R, sea ϕu:R2−→ R una función de densidad de probabilidad con media un vector u∈R2dado y sea T:R2−→ R2 (x, y)7−→ (x, y −x). Con todos estos ingredientes, podemos definir la superficie de persistencia de B como la función: ρB:R2−→ R2 z7−→ ρB(z) = X u∈T(B) g(u)ϕu(z). Dado el conjunto de puntos de un diagrama de persistencia B(sin los puntos que tienen segunda coordenada infinita), aplicamos la transformación lineal T para pasar de coordenadas nacimiento-muerte en Ba coordenadas nacimientopersistencia en T(B). Después, obtenemos la superficie de persistencia como la suma ponderada de funciones de densidad centradas en los puntos. La superficie de persistencia ρBse reduce a un vector finito, fijando una rejilla a lo largo de una región de interés acotada del plano (por ejemplo, la región más pequeña que contiene todos los puntos de T(B)) y obtenemos una imagen de persistencia integrando ρBen todos los cubos de la rejilla. Definición 2.18. Sea Bel conjunto formado por todos los puntos de un diagrama de persistencia que tienen multiplicidad positiva y no tienen persistencia infinita y ρBuna imagen de persistencia del diagrama. Sea {I1, I2, . . . , In}unos píxeles que cubren como rejilla una región de interés acotada del plano. Así la imagen de persistencia de Bes el vector I(ρB)I1=Z ZI1 ρBdydx, I(ρB)I2=Z ZIn ρBdydx, . . . , I(ρB)In=Z ZIn ρBdydx. 2. conceptos previos 25 Figura 2.7: Podemos ver el proceso de obtención de una imagen de persistencia a partir de una imagen en escala de grises. En la imagen de la izquierda vemos una imagen en escala de grises como datos de input, a su izquierda vemos un diagrama de persistencia Basociado de una dimensión k∈ {0,1}(sin los puntos que tienen persistencia infinita) y a la izquierda de Bvemos T(B), i.e., el diagrama Bcambiando de coordenadas nacimiento-muerte a nacimiento-persistencia. La cuarta imagen por la izquierda es la superficie de persistencia asociada a B. En la quinta imagen por la izquierda podemos ver varias imágenes de persistencia diferenciadas por la precisión que se usa para integrar la superficie de persistencia anterior. Además, en la quinta foto por la izquierda podemos ver una imagen de persistencia aplanada, formando un vector. Imagen sacada de [15]. En la figura 2.7 podemos ver de manera gráfica la obtención de una imagen de persistencia a partir de una imagen en escala de grises. Las imágenes de persistencia dependen de 3parámetros: Resolución de la imagen: Mide el tamaño de los píxeles considerados al integrar la superficie de persistencia. Cuánto más resolución, menos tamaño de los píxeles y más dimensiones va a tener la imagen de persistencia. Simulaciones en [15], muestran que en métodos de clasificación a partir de imágenes de persistencia el porcentaje de acierto parece robusto frente la elección de resolución. Esto indica que, aparentemente, reducir la dimensionalidad de las imágenes de persistencia al reducir la resolución no afecta mucho a la calidad de la información extraída por las imágenes de persistencia. La distribución ϕu. Normalmente, se suele emplear como ϕuuna distri- 26 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos bución normal centrada en u= (ux, uy)∈T(B)y con matriz de varianzascovarianzas con σ2en los elementos de la diagonal y 0en el resto, i.e., se suele considerar ϕu=1 2πσ2e−[(x−ux)2+(y−uy)2]/2σ2 para un u∈T(B). Si se escoge la distribución normal, queda por escoger otro parámetro, σ, la desviación típica de cada variable. Experimentos en [15], muestran que en métodos de clasificación a partir de imágenes de persistencia el porcentaje de acierto es poco sensible frente la elección de σ. Función de ponderación g. Para asegurar la estabilidad (ver 5.1), se exige que gsea continua, diferenciable a trozos y que g(x, 0) = 0 para todo x∈R. Estas funciones, además, buscan ponderar más a los puntos con más persistencia, los más alejados de la diagonal en los diagramas, ya que son los que más información topológica sobre la estructura de los datos nos aporta (y es menos probable que sea ruido topológico). Una opción simple es que gdependa solo de la segunda coordenada (la de persistencia) y que, para ponderar los puntos del diagrama más cuanta más persistencia tengan se puede considerar una función no decreciente respecto de la segunda coordenada. En connivencia con lo anterior se puede usar g(x, y) = yoguna función lineal a trozos en la variable yque valga 0si y= 0 y1para el punto del diagrama con más persistencia. 2.3 Modelos de aprendizaje automático. Clasificación supervisada Un modelo de aprendizaje automático es un modelo en que una máquina (un ordenador normalmente) obtiene un conocimiento a partir de unos datos que se le aportan y el desempeño de la máquina mejora cuanto más datos se le aportan. "En el aprendizaje automático, un computador observa datos, construye un modelo basado en esos datos y utiliza ese modelo a la vez como una hipótesis acerca del mundo y una pieza de software que puede resolver problemas" ( [16, p. 651]). Podemos clasificar el aprendizaje automático según los datos "a priori" que tenemos. El aprendizaje supervisado es un tipo de aprendizaje automático en que en 2. conceptos previos 27 los datos existen una o más variables a predecir (variables dependientes) en función de un conjunto de variables predictoras (variables independientes). Por tanto una manera de evaluar la calidad de los modelos de aprendizaje supervisado es comparar los valores reales de las variables a predecir con los valores de las predicciones, de manera que cuanto más se parezcan mejor será el modelo. El aprendizaje no supervisado es un tipo de aprendizaje automático en que no existen un conjunto de variables a predecir y la máquina se dedica a aprender patrones que tengan los datos. Por tanto, como esos patrones de los datos no se saben "a priori", ya que no son variables de los datos, no se puede evaluar la calidad de los modelos comparándolos con unos valores "reales", sino que la evaluación de la calidad de los modelos se centra en medir cómo de bien los patrones obtenidos explican los datos. Los métodos de aprendizaje automático supervisado se pueden dividir a su vez según el tipo de variable dependiente: Regresión: La variable a predecir es cuantitativa, i.e., numérica (tanto continua como discreta). Clasificación: La variable a predecir es cualitativa y por tanto no se pueden hacer operaciones matemáticas con ella. En algún caso las variables cualitativas se sustituyen por el valor de un número (etiquetas numéricas), pero al representar categorías (raza, religión, estado civil, ...) no se pueden hacer operaciones con ellas. En caso de que la variable a predecir sea dicotómica (tenga solo dos valores o categorías posibles) estamos ante clasificación dicotómica o binaria. En este trabajo emplearemos métodos supervisados de clasificación dicotómica. En la segunda parte del método descrito en el capítulo 3 se emplea una máquina vector soporte, que explicaremos en la sección 2.4. En la segunda parte del método descrito en el capítulo 4, que describiremos en la sección 2.5. En cada método se ha usado uno de esos métodos de aprendizaje automático porque eran los mejores (según el procedimiento de la sección 2.3.2) de un conjunto de métodos candidatos. Antes de ello vamos a explicar cuestiones que son comunes a los métodos de clasificación dicotómicos. 28 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos 2.3.1 Sobreajuste, división de conjuntos, validación cruzada. Cuando trabajamos con un conjunto de datos, ese conjunto de datos suele ser una muestra de una población (por ejemplo: los resultados de una encuesta es una muestra respecto de la población total a la que se dirigía la encuesta) y por tanto la información que saquemos de los datos vamos a tomarla como una hipótesis de toda la población. Por lo tanto, queremos que los modelos que obtengamos de los datos también funcionen en la población. Aquí es donde entra el sobreajuste, que significa que el modelo funciona muy bien en el conjunto de datos con el que la máquina aprende, pero no tan bien en el total de la población. Para evitar que el modelo que seleccionemos (en este caso de clasificación dicotómica) se ajuste muy bien a los datos pero no se pueda generalizar, dividimos el conjunto de datos que trabajamos en dos de manera aleatoria: Conjunto de entrenamiento: Suele contener entre el 70 y el 80 por ciento de los datos. Sobre este conjunto se construyen los modelos. Conjunto test: Suele contener entre el 20 y el 30 por ciento de los datos. Sobre este conjunto se evalúan los modelos construidos sobre el conjunto de entrenamiento. Así, con esta división no solo buscamos un modelo que se ajuste bien a los datos con los que se contruye, sino que se pueda generalizar. La capacidad de generalización la medimos en el conjunto test, ya que ese conjunto no se ha usado para construir el modelo. Cuando construimos un modelo a partir de un método de aprendizaje automático, ese modelo no solo va a depender de los datos del conjunto de entrenamiento sino que puede depender de uno o más hiperparámetros, como por ejemplo en una regresión polinomial sería el grado del polinomio del modelo o en una regresión múltiple cuantas variables tiene el modelo. Para escoger estos hiperparámetros tenemos un problema parecido al que planteábamos antes, saber a partir del conjunto de entrenamiento cuáles serán los mejores hiperparámetros en el conjunto test. Además muchas veces con unos valores para los hiperparámetros obtenemos un modelo muy complejo que se sobreajusta al conjunto de entrenamiento (se ajusta muy bien al conjunto de entrenamiento pero tiene peor rendimiento en el conjunto test) y con otros valores obtenemos modelos muy simples que no tienen buen rendimiento ni en el conjunto de entrenamiento ni en el test. 2. conceptos previos 29 Un método que se emplea para escoger los hiperparámetros correctos es la validación cruzada. Con este método se divide el conjunto de entrenamiento Cen kpartes C1, . . . , Ckde igual tamaño. Nos quedamos con unos posibles valores para los hiperparámetros a escoger. Para cada parte idel conjunto de entrenamiento, i∈ {1, . . . , k}construimos sobre C\Ciun modelo para cada valor posible de los hiperparámetros y evaluamos el modelo en el conjunto de datos Ci. Finalmente, nos quedaremos con los valores para los hiperparámetros que mejor evaluación obtenga (considerando la media de las evaluaciones en los conjuntos C1, . . . , Ck. Normalmente, se suelen considerar k= 5 o10. Cómo evaluar un modelo de clasificación binaria se verá en la subsección 2.3.3. 2.3.2 Elección del mejor método de aprendizaje automático con validación cruzada anidada LTO Supongamos que tenemos varios métodos de aprendizaje automático, y queremos escoger que método funciona mejor con un tipo de datos que tenemos. Para ello, vamos a realizar un tipo de validación cruzada algo distinta al de la sección 2.3.1 en el que vamos a testear los métodos sobre todo el conjunto. El método que vamos a usar, propuesto por Thi Lan Anh Dinh et al. en [17], es la validación cruzada anidada dejando dos fuera (leave-two-out, LTO). Primero, dividimos el conjunto de datos en ksubconjuntos del mismo tamaño: C1, . . . , Ck. Después, para cada iteración i∈ {1, . . . , k −1}, construimos un modelo en el conjunto de entrenamiento C1∪ · · · ∪ Ci−1∪Ci+1 ∪ · · · ∪ Ckcon los hiperparámetros que mejor evaluación obtenga en el conjunto de validación Ci+1 y evaluamos ese modelo en el conjunto test Ci. En la iteración i=khacemos lo mismo solo que con el conjunto de entrenamiento C2∪ · · · ∪ Ck−1, de validación Cky test C1. Nos quedamos con el método que obtenga la mejor evaluación media. Observar que cada modelo que construimos y evaluamos en cada iteración se construye usando hiperparámetros fijados únicamente para esa iteración. Se suele considerar k= 5 o10. Una vez hayamos escogido el mejor método para nuestros datos, podemos considerar, para dicho método, la validación cruzada (VC) explicada en la subsección 2.3.1, i.e., dividimos en conjunto test y entrenamiento una única vez y fijamos los hiperparámetros solo una vez por VC dividiendo el conjunto de entrenamiento. Así, obtendríamos el método entrenado con unos hiperparámetros fijos. 30 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Cuadro 2.1: Una matriz de confusión para la clasificación ˆyde y. Valores Reales y= 1 y=−1 Valores ˆy= 1 Verdaderos Pos. (VP) Falsos Pos. (FP) VP+FP Predichos ˆy=−1Falsos Neg. (FN) Verdaderos Neg. (VN) FN+VN VP+FN FP+VN n=VP+FP+FN+VN 2.3.3 Evaluación de la calidad de una clasificación supervisada binaria Sea un conjunto de datos (x1, y1),...,(xn, yn), con xi∈Rpeyi∈ {−1,+1} ∀i∈ {1, . . . , n}donde las primeras pcoordenadas son las variables predictoras y la última coordenada es la variable dicotómica ya predecir. Supongamos que hemos aplicado un método de aprendizaje automático en el que hemos obtenido una predicción ˆyde los valores de la variable y, denotando para cada i∈ {1, . . . , n} la predicción de yicomo ˆyi, entonces un primer método de evaluar la calidad de la clasificación ˆyde yes la realización de una tabla de confusión. La matriz de confusión consiste en una matriz que cuenta el número de veces que ˆytoma cada valor de {−1,+1}cuando yvale un valor {−1,+1}. Es decir, una matriz con la frecuencia absoluta conjunta de yeˆy. En el cuadro 2.1 podemos ver como son las matrices de confusión y lo que significa cada elemento de la matriz. Así, una clasificación binaria será mejor cuanto los elementos que tiene en la diagonal sean más grandes y más pequeños los elementos fuera de la diagonal. Definición 2.19. A partir de la tabla de confusión se pueden calcular la siguientes medidas de la bondad de la clasificación: Exactitud (o accuracy en ingés): Proporción de datos clasificados correctamente, i.e., la suma de la diagonal de la matriz de confusión entre la suma del total de elementos de la matriz: AC := V P +V N n. Sensibilidad: Proporción de positivos (en este caso de datos con y= 1) clasificados correctamente, i.e., SE := V P V P +FN . 2. conceptos previos 31 Cuadro 2.2: Matriz de confusión del ejemplo 2.4. Valores reales y= 1 y=−1 Valores ˆy= 1 90 10 100 predichos ˆy=−10 0 0 90 10 100 Precisión: Proporción de datos positivos entre el total de los clasificados como positivos, i.e., PR := V P V P +FP . Especificidad: Proporción de negativos (en este caso de datos con y=−1) clasificados correctamente, i.e., ES := V N FP +V N . Todas estas medidas valen siempre entre 0y1y el modelo será mejor cuando estas medidas más se acerquen a 1. La medida o medidas que son más apropiadas para valorar el modelo depende del objetivo que tengamos con el modelo (clasificar bien los datos en general, clasificar bien los positivos/negativos, tener pocos falsos positivos, ...) Puede ser que dos medidas nos den un resultado muy bueno pero, sin embargo, otra medida tenga un resultado muy pobre. Esto puede ocurrir cuándo el número de elementos de una clase y otro son muy distintos. Lo podemos ver en este simple ejemplo. Ejemplo 2.4.Un método de clasificación que tenga la matriz de confusión del cuadro 2.2 tiene: AC = 0.9,SE = 1 yES = 0. Que se clasifiquen tan mal los negativos no afecta mucho a la proporción total de datos clasificados correctamente (exactitud), ya que solo el 10 por ciento de los datos son negativos. 32 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos 2.4 Máquina Vector Soporte Una máquina vector soporte (MVS) es un método de aprendizaje automático que clasifica una variable de un conjunto de datos (p+ 1)-dimensionales según si el vector p-dimensional de sus variables predictoras estén a un lado u otro de un hiperplano (óptimo). Esto significa que, dado un conjunto de datos (x1, y1),...,(xn, yn), con xi∈Rpeyi∈ {−1,+1} ∀i∈ {1, . . . , n}, la MVS construye una forma de predecir si yi= +1 o−1, en función de los valores de su vector xiy, para ello, separaremos los vectores xien dos grupos separados por un hiperplano en Rp. La máquina vector soporte se define inicialmente para clasificar una variable dicotómica, pero se puede generalizar a variables cuantitativas. En este trabajo consideraremos solo variables dicotómicas donde yi= 1 si ipertenece a la clase C1oyi=−1si ipertenece a la clase C2. Primero, tratamos el caso de que los datos sean separables linealmente en Rpsegún su clase, luego, generalizaremos a datos casi linealmente separables introduciendo términos de error, y finalmente, generalizaremos el modelo para poder hacer separaciones no lineales proyectando los datos a dimensiones superiores. 2.4.1 Datos separables linealmente Sean (x1, y1),...,(xn, yn), con xi∈Rpeyi∈ {−1,+1} ∀i∈ {1, . . . , n}. Si suponemos que los datos se pueden separar linealmente, entonces existe al menos un hiperplano xtβ+β0= 0 con β∈Rp, β0∈Rtal que para todo i∈ {1, . . . , n} se cumple que o bien xt iβ+β0≥0si y= 1 o bien xt iβ+β0≤0eyi=−1. Ahora, vamos a intentar que los datos estén lo más separados del hiperplano separador. Dado que la distancia de un dato xial hiperplano xtβ+β0= 0 es: |xt iβ+β0| ||β|| ; si y= 1, como xt iβ+β0≥0, vamos a intentar maximizar xt iβ+β0 ||β|| y si y=−1, como xt iβ+β0≤0, vamos a intentar minimizar xt iβ+β0 ||β|| . Entonces obtenemos el siguiente problema de optimización: 2. conceptos previos 33 m´ax ||β||=1M s.a. yi(xt iβ+β0)≥M∀i∈ {1, . . . , n} (2.1) Con este problema estamos intentando maximizar la distancia del dato más cercano al hiperplano. Como los vectores normalizados en Rpincluyen todas las direcciones que hay en Rppodemos restringirnos al caso ||β|| = 1 sin pérdida de generalidad. Si consideramos todos los vectores β∈Rpsin normalizar modificando los β0 por ||β||β0las restricciones del problema son equivalentes a yi ||β||(xt iβ+β0)≥M∀i∈ {1, . . . , n}. Estas condiciones son equivalentes a yi(xt iβ+β0)≥M||β||. Como dados unos β∈Rp, β0∈Rque satisfagan la última desigualdad ∀i∈ {1, . . . , n}, si los multiplicamos por un número real positivo cualquiera se siguen satisfaciendo las desigualdades y por tanto podemos tomar arbitrariamente ||β|| =1 M. Así, reescalando otra vez el β0, aplicando sobre la función objetivo la función monótona decreciente f(x) = 1 2x2y cambiando maximizar por minimizar, (2.1) es equivalente a: m´ın 1 2||β||2 s.a. yi(xt iβ+β0)≥1∀i∈ {1, . . . , n} (2.2) Con las restricciones estamos creando un par de rectas paralelas xtβ+β0= 1 yxtβ+β0=−1donde si yi= 1 entonces xt iβ+β0≥1yyi= 1 entonces xt iβ+β0≤ −1para todo i. Es decir, estamos creando dos rectas paralelas donde no haya ningún dato y donde al otro lado de cada recta solo haya puntos de una clase. La distancia entre esas rectas paralelas, el ancho de ese margen entre clases, es 2 ||β||, que es lo que estamos maximizando al minimizar ||β||2. 2.4.2 Datos cuasi-separables linealmente En el caso de que los conjuntos no sean separables linealmente, vamos a introducir al problema (2.2) un término de error para cada punto i∈ {1, . . . , n}: 3 Estado del arte En este capítulo describiremos el método de Muszynski, G. et al. ([7]) para detectar si existe o no un río atmosférico en una región de la Tierra en un determinado momento de tiempo. Este método extrae patrones topológicos sobre los datos, en particular, la evolución de las componentes conexas del subgrafo de nivel superior. Luego, con aprendizaje automático, se intenta diferenciar qué patrones se corresponden con ser RA, dando resultado a una clasificación de las imágenes como RA o como no RA. Para el aprendizaje automático se emplean unas etiquetas provenientes de otro método, que también se usan para medir la calidad de la clasificación. En la sección 3.1 se describen con qué conjuntos de datos se ha trabajado y en la sección 3.2 se muestra un esquema del método. A continuación, detallamos la primera etapa del método, donde se usa ATD, en la sección 3.3 y detallamos la segunda etapa, la del aprendizaje automático, en la sección 3.4. Finalmente, se muestran y analizan los resultados que se obtienen en la sección 3.5. 3.1 Datos empleados Los conjuntos de datos que se emplean vienen tanto de la simulación de modelos climáticos generada por la versión 5.1 de la Community Atmosphere Model (CAM5.1) ([20]), como del producto del reanálisis de la segunda versión del Modern-Era Retrospective Analysis for Research and Applications (MERRA-2) ([21]). El reanálisis es un proceso a través del cual un método fijo de asimilación de datos (que depende de un modelo meteorológico de predicción) se usa para adquirir reprocesamientos consistentes de observaciones meteorológicas. De estos datos se emplea la columna que nos indica el DSVA. 42 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Cuadro 3.1: Tabla con los datos que se emplean, tanto provenientes de simulaciones como de reanálisis, y con sus períodos y resoluciones. Modelo climático Período Resolución temporal Resolución espacial CAM5.1 1979 −2005 3h y 24h 25 km CAM5.1 1979 −2005 3h y 24h 100 km CAM5.1 1979 −2005 3h y 24h 200 km MERRA-2 1980 −2017 3h 50 km Aunque TDLVA es otra variable también usada en la identificación de ríos atmosféricos, no es observable por satélite, al contrario que la DSVA. Sin embargo, aunque no se contempla en esta memoria, también se pueden hacer modelos de identificación de RA basados en TDLVA. Estos modelos dan imágenes en escala de grises de dimensión 2, que son observaciones en un instante de tiempo del valor de DSVA a lo largo de una rejilla que tiene puntos a lo largo de un mapa de una región. Estas imágenes en escala de grises vienen representadas por matrices. En este caso se emplea como región el área del pacífico que va desde las islas de Hawái en EE. UU. hasta la costa oeste de los EE. UU. (la región que se ve en la figura 1.1). Se utilizan datos con distintas precisiones temporales (distancia de tiempo entre los instantes de tiempo en los que observamos la variable) y espaciales (distancia entre los puntos de la rejilla contiguos). Los modelos de datos empleados, precisiones y períodos de tiempo correspondientes se pueden ver en el cuadro 3.1. Para que se vea el formato y el tamaño de los datos, veamos una submatriz 7×7de una matriz (73 ×112) de uno de los conjuntos usados, el proveniente del MERRA-2 con precisión espacial de 50km y temporal de 3h:            23.658932 23.036139 22.123043 21.430202 22.051882 23.914698 25.141624 21.240788 20.877882 21.092453 22.048023 23.707134 25.040762 25.717901 20.433863 21.380583 22.717028 24.17173 25.391476 25.826097 25.939281 22.098997 23.790562 25.248415 25.969824 26.210562 26.03678 25.719465 24.092237 25.286072 26.315365 26.481539 26.273636 26.1369 26.369896 24.930016 25.768875 26.599392 26.673508 26.62879 27.229403 28.029497 25.454836 25.905527 26.493378 27.133888 28.030231 28.798422 28.613619            3. estado del arte 43 Por otra parte, como se tiene que entrenar un clasificador de aprendizaje automático no supervisado, se necesita que las observaciones estén etiquetadas como RA o no RA, es decir, sobre si en cada imagen en escala de grises hay un río atmosférico o no. Las etiquetas se obtienen aplicando a los datos el modelo Toolkit for Extreme Climate events Analysis (TECA) ([22]). Este método clasifica según un criterio fijo basado en un umbral. Se asume que las etiquetas que aporta el TECA como "etiquetas verdaderas". Los datos que se usan están disponibles en [23]. 3.2 Esquema del método Figura 3.1: Esquema del método. Imagen sacada de [7]. El método consta de dos etapas: Etapa 1: Se extraen características topológicas de las imágenes con Análisis Topológico de Datos. Se emplea el algoritmo Union-Find para sacar información sobre la evolución de las componentes conexas en unos conjuntos de nivel. Estas características topológicas serán el input del clasificador de la segunda etapa. Etapa 2: Clasificamos las características previas con un método aprendizaje automático, la Máquina Vector Soporte (MVS). La clasificación se divide en dos pasos: Primero, se entrena el modelo para separar las características topológicas en dos clases (RA y no RA) en el conjunto de entrenamiento y después se le aplica el modelo al conjunto test. Posteriormente, se evalúa 44 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos la clasificación en términos de exactitud, sensibilidad y precisión. Como estamos realizando aprendizaje supervisado, tanto para entrenar el modelo como para evaluar sus resultados, vamos a emplear las etiquetas que nos aporta el modelo TECA. El método tiene los siguientes datos de entrada y de salida. Datos de entrada: nimágenes en escalas de grises bidimensionales con datos en ninstantes de tiempo sobre la densidad superficial de vapor de agua (DSVA) en kgm−2. Cada imagen en escala de grises viene representada por una matriz. Para cada imagen tenemos una etiqueta (las obtenidas del modelo TECA) que nos indica a que clase pertenece (RA o no RA) esa observación. Datos de salida. Etiquetas binarias. Si es RA,RA = 1, si no lo es, RA = 0. 3.3 Primera etapa, análisis topológico de datos En esta etapa solo se emplean las imágenes sin las etiquetas (es decir, las coordenadas y el valor de DSVA). Así, los datos se pueden representar como una función definida en una rejilla, determinada por las coordenadas de los datos, al conjunto de los reales, donde la imagen de las coordenadas es el valor del DSVA en ese punto. Como el valor de DSVA es positivo y nunca supera su máximo Len la rejilla (existe el máximo en la rejilla al ser un conjunto finito de puntos) entonces esa función tiene su imagen dentro de [0, L], donde en este caso L= 60kgm−2 Esta función sería la restricción de una función DSVA: [a, b]×[c, d]−→ [0, L′] definida en toda la región del estudio (y no únicamente en la rejilla) donde a,b, cydrepresentarían los límites izquierdo/derecho/superior/inferior de la región y donde la función asigna a cada punto su valor de la densidad superficial de vapor de agua. L′sería el valor máximo de DSVA en la región. Por otra parte, podemos modelar la rejilla como un grafo donde cada punto está unido por aristas a los puntos que están inmediatamente a su izquierda, 3. estado del arte 45 Figura 3.2: Representación de la rejilla como un grafo donde cada punto tendría como aristas los vecinos inmediatamente superiores, inferiores, a su izquierda y derecha. Imagen sacada de [7]. derecha, arriba y abajo. Así, cada punto tendría 4 vecinos (salvo los puntos de la frontera). Podemos ver una representación del grafo en la figura 3.2. Definimos el conjunto de nivel superior de f(en la región [a, b]×[c, d]) como los puntos de la región que tienen DSVA mayor que un valor t: (DSVA)−1[t, +∞) = {(x, y)∈[a, b]×[c, d]|f(x, y)≥t}. Interesa saber como evoluciona la forma de las zonas de la región con DSVA mayor que cierto umbral a medida que va disminuyendo el umbral, ya que los ríos atmosféricos tienen una forma característica (son bandas alargadas y estrechas dónde hay un alto valor de DSVA). En particular, se va a ver como evolucionan las componentes conexas de (DSVA)−1[t, +∞). Como solo se saben los valores en la rejilla, se va a ver como evoluciona la conexidad de los subconjuntos de nivel superior de DSVA restringida a la rejilla. Es decir, son subgrafos dónde se considera que: Un nodo está en el subconjunto de nivel si su valor de DSVA es mayor o igual que el umbral. 46 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Una arista está en un subconjunto de nivel superior si lo están los dos puntos (nodos) de la rejilla que une. En el grafo se considera que dos puntos están en la misma componente conexa si existe un camino que los une. Se va a calcular el conjunto de nivel superior de fen la rejilla para los valores t∈ {1,1 + p, 1+2p, . . . , L = 60 kgm−2}para el umbral (siendo p > 0un paso entre umbrales considerados) y se va a medir cómo evolucionan las componentes conexas del conjunto según vayamos reduciendo el valor del umbral. Cuando t > L este conjunto es vacío, cuando tse va reduciendo aparecen componentes conexas que a medida que se sigue reduciendo tse van a sumergir unas en otras. Solo interesan las componentes conexas que intersequen ciertas subregiones del mapa. En el caso particular de este estudio, como se está tratando en el pacífico noreste, solamente interesan las componentes conexas que intersequen la latitud de Hawái y el continente americano. En la figura 3.3 se pueden ver las subregiones que interesan (la costa americana rodeada de verde y las latitudes de Hawái en amarillo). Si alguna componente conexa no interseca la subregión de interés, se actúa como si no existiese esa componente conexa, ya que esa componente conexa no puede forma un patrón de río atmosférico, ya que no va desde los trópicos hasta el continente. 3.3.1 Gráficos de evolución Para obtener las características topológicas que luego se usan como input en la máquina vector soporte se emplea un algoritmo basado en el algoritmo Union- Find data structure ([24]). Con este algoritmo se obtiene el número de puntos de la rejilla que están en una componente conexa a medida que el umbral disminuya (se toman como valores para el umbral 1,1 + p, 1+2p, . . . , L = 60 kgm−2,p > 0)1. Para determinar en qué componente conexa se cuentan los puntos de la rejilla se emplea la Elder Rule. La Elder Rule consiste en que cuando dos componentes conexas se unen, se considera que la que ha aparecido antes (la que apareció 1En la implementación hecha por Muszynski, G. et al. se ha considerado un paso p= 1. 3. estado del arte 47 Figura 3.3: Una ilustración de la búsqueda de regiones conectadas de RA en unas subregiones específicas de interés. En este ejemplo, las subregiones que interesan son la costa americana rodeada de verde y las latitudes de Hawái en amarillo. (a) El cuadrado rojo indica el lugar donde dos componentes conexas del conjunto de nivel superior no están conectadas para un valor del umbral (DSV A =t1). (b) Con nuevo umbral (DSV A =t2), donde t2< t1, las dos componentes conexas se unen haciendo que el conjunto de nivel superior tenga forma de RA. El DSVA (en kgm−2) mostrado en este ejemplo proviene del modelo CAM5.1. Imagen sacada de [7]. con un umbral más grande) es la que absorbe a la otra y sobrevive. Por tanto, con la Elder Rule estamos midiendo el número de elementos que tiene para cada umbral la primera componente que aparece (ya que es la única que cuando dos componentes conexas se unen siempre sobrevive). El algoritmo descarta las componentes conexas que no interseca las subregiones de interés y no las emplea en sus cálculos. Esta información que nos aporta el algoritmo la podemos representar en gráficos de evolución. Los gráficos muestran el número de puntos rejilla en una componente conexa (la que aparece antes) según el valor de DSVA. Estos gráficos muestran el número normalizado de puntos de la rejilla, es decir, el porcentaje de puntos de la rejilla, que están en una componente conexa para cada valor umbral. Se normalizan los puntos para así poder comparar gráficos de rejillas con distinta resolución espacial (por tanto, distinto número total de puntos en la rejilla). En la figura 3.4 podemos ver tres gráficos de evolución para el mismo conjunto de datos: a la izquierda podemos ver los gráficos de evolución de 100 observaciones etiquetadas como ríos atmosféricos y la media de todas las observaciones 48 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos etiquetadas como RA; en el centro los gráficos de evolución de 100 observaciones etiquetadas como no RA y la media de todas las observaciones etiquetadas como no RA; y en el centro se comparan las dos medias previas. Estos gráficos medios son muy parecidos, casi indiferenciables para un ojo humano, pero son distintos, y es algo de lo que se encargará después la máquina vector soporte. Figura 3.4: Gráficos de evolución normalizados de 100 observaciones y la media (rojo) de todas las observaciones. En la parte de arriba, a la izquierda, vemos el grafo de evolución para observaciones etiquetadas por TECA como RA (azul) y, en la derecha, para etiquetadas como no RA. En la parte de abajo, comparamos ambas medias. Las observaciones se corresponden con los datos sacados del modelo de simulación CAM 5.1 con precisión espacial de 200km y observaciones diarias. Imagen sacada de [7]. De estos gráficos de evolución se puede extraer por cada observación un vector kdimensional donde la coordenada i-ésima del vector se corresponde con el número normalizado de puntos en una componente conexa considerando el umbral i-ésimo. Así, se obtendría, juntando todas las observaciones, una matriz con n filas y kcolumnas, donde nes el número de observaciones (instantes de tiempo en el que se tiene una imagen de la variable DSVA en el mapa) y kes el número de umbrales que se consideran (por tanto, el tamaño de las características 3. estado del arte 49 topológicas extraídas de cada observación). Esta matriz es el input de la máquina vector soporte descrito en la siguiente sección. En este caso, como tenemos que L= 60kgm−2es el valor máximo de DSVA en la rejilla, se consideran los umbrales 1,2,...,60 = L, y por tanto, se considera k= 60 y un paso entre umbrales de p= 1. Si redujésemos el paso, aumentaríamos la dimensión de los datos. En la figura 3.5 podemos ver un esquema de cómo se genera la matriz input. Figura 3.5: Esquema de cómo se crea la matriz input para la MVS. Imagen sacada de [7]. 3.3.2 Ejemplo simple del funcionamiento del algoritmo En esta sección vamos a mostrar un ejemplo simple de como funciona el algoritmo. Vamos a suponer que nuestra matriz de datos (que nos da el valor de DSVA a lo largo de la rejilla) es una matriz 5×5. Modelamos entonces la rejilla como un grafo con 25 nodos que están unidos por aristas con los nodos justamente inferiores, superiores, a su izquierda y a su derecha. Además, por simplificar el 56 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos En resumen, el modelo de detección de RA obtiene una alta exactitud (0.77- 0.91) aplicándolo a diversos conjuntos de datos con distinto origen, precisión temporal y precisión espacial. Esto ilustra que la combinación de análisis topológico de datos y aprendizaje automático es una estrategia efectiva y eficiente para detectar RA en conjuntos de datos climáticos grandes sin tener que fijar previamente umbrales subjetivos. Además, el modelo de aprendizaje automático está sesgado por las etiquetas empleadas provenientes de TECA, ya que hemos asumido que estas etiquetas son "verdad" y hemos entrenado y calificado nuestro modelo con ellas. Queda pendiente de comprobar y usar el método con distintos etiquetados provenientes de otros métodos. 4 Método propuesto En este capítulo proponemos un método que mejora el de Muszynski, G. et al. [7]. En el capítulo 5 compararemos los resultados y veremos que el método que proponemos mejora en términos de estabilidad y de resultados numéricos. Al igual que el método de Muszynski, G. et al. con este método vemos si existe o no un río atmosférico en una región de la Tierra en un determinado momento de tiempo. Este método extrae patrones topológicos sobre los datos, en particular, la evolución de las componentes conexas y de los huecos en los complejos cúbicos filtrados asociados a las imágenes en escala de grises. Luego, con aprendizaje automático, se intenta diferenciar qué patrones se corresponden con ser RA, dando resultado a una clasificación de las imágenes como RA o como no RA. Para el aprendizaje automático se usan unas etiquetas provenientes de otro método, que también se usan para medir la calidad de la clasificación. La principal diferencia de nuestro método con el de Muszynski, G. et al. es que no solo extraemos información sobre la evolución de las componentes conexas sino que también sobre la evolución de los huecos. Es por ello que esperamos que mejore los resultados de Muszynski, G. et al.. Además: las etiquetas que usamos en el aprendizaje automático (provenientes de TECA-BARD) son distintas a las que usan Muszynski, G. et al. en su método (provenientes de TECA) y empleamos otro método de aprendizaje automático. En la sección 4.1 se describen con qué conjuntos de datos se ha trabajado y en la sección 4.2 mostramos un esquema del método propuesto y qué diferencias y similitudes tiene con el anterior. A continuación, detallamos la primera etapa de nuestro método, donde se usa ATD, en la sección 4.3 y detallamos la segunda etapa, la del aprendizaje automático, en la sección 4.4. Finalmente, hablamos sobre el código y la implementación a ordenador en 4.5 y los resultados que se obtienen se muestran y analizan en 4.6. 58 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Cuadro 4.1: Tabla con los datos que se emplean, provenientes de simulaciones como de reanálisis, y con sus períodos y resoluciones. Modelo climático Período Resolución temporal Resolución espacial MERRA-2 1980 −2017 3h y 24 h 50 km 4.1 Datos empleados Los conjuntos de datos que se emplean vienen del producto del reanálisis de la segunda versión del Modern-Era Retrospective Analysis for Research and Applications (MERRA-2) ([21]). El reanálisis es un proceso a través del cual un método fijo de asimilación de datos (que depende de un modelo meteorológico de predicción) se usa para adquirir reprocesamientos consistentes de observaciones meteorológicas. Al igual que en el capítulo anterior, de estos datos se emplea la columna que nos indica el DSVA y no se utiliza TDLVA, que es otra variable que se emplea en la detección de ríos atmosféricos, ya que no es observable por satélite, al contrario que la DSVA. Este modelo da imágenes en escala de grises de dimensión 2, que son observaciones en un instante de tiempo del valor de DSVA a lo largo de una rejilla que tiene puntos a lo largo de un mapa de una región. Estas imágenes en escala de grises vienen representadas por matrices. Al igual que en el método descrito en el capítulo 4 se usa como región el área del pacífico que va desde las islas de Hawái en EE. UU. hasta la costa oeste de los EE. UU. (la región que se ve en la figura 1.1). Se emplean datos con distintas precisiones temporales (distancia de tiempo entre los instantes de tiempo en los que observamos la variable). Los modelos de datos empleados, precisiones y periodos de tiempo correspondientes se pueden ver en el cuadro 4.1. Para que se vea el formato y el tamaño de los datos provenientes del MERRA- 2, veamos una submatriz 7×7de una matriz (73×112) de los conjuntos utilizados: 4. método propuesto 59            23.658932 23.036139 22.123043 21.430202 22.051882 23.914698 25.141624 21.240788 20.877882 21.092453 22.048023 23.707134 25.040762 25.717901 20.433863 21.380583 22.717028 24.17173 25.391476 25.826097 25.939281 22.098997 23.790562 25.248415 25.969824 26.210562 26.03678 25.719465 24.092237 25.286072 26.315365 26.481539 26.273636 26.1369 26.369896 24.930016 25.768875 26.599392 26.673508 26.62879 27.229403 28.029497 25.454836 25.905527 26.493378 27.133888 28.030231 28.798422 28.613619            Por otra parte, como se tiene que entrenar un clasificador de aprendizaje automático no supervisado, se necesita que las observaciones estén etiquetadas como RA o no RA, es decir, sobre si en cada imagen hay un río atmosférico o no. Las etiquetas se obtienen aplicando a los datos el modelo TECA Bayesian Atmospheric River Detector (TECA-BARD) ([25]). Este método automático de detección de RA se ha construido a partir de datos etiquetados manualmente por expertos. Se asume que las etiquetas que aporta el TECA-BARD como "etiquetas verdaderas". En comparación, vamos a usar un conjunto de datos de los que usan Muszynski, G. et al. en su método ([7]) y además vamos a usar ese conjunto de datos, pero con solo un dato cada día. Sin embargo, no se van a emplear las mismas etiquetas para ese conjunto de datos, sino que las etiquetas que vamos a utilizar son provenientes de un método distinto. Esta diferencia simplemente es debida a las dificultades de acceso a los datos. A pesar de ello, los resultados son comparables pero teniendo en cuenta ese matiz. Los datos que se utilizan están disponibles en [23]. Las etiquetas para los datos empleados 1están disponibles en el siguiente link: https://www.earthsystemgrid. org/dataset/ucar.cgd.ccsm4.artmip.tier1.cascade_bard_v1.atm.other.3hourly_ inst.html 1Se han procesado los datos de una rejilla a lo largo de la región (por cada instante de tiempo) con 1en los puntos dónde había un río atmosférico y 0si no. Hemos obtenido, a partir de cada matriz (rejilla) un único número por cada rejilla que nos indicase si en la región de estudio había un río atmosférico que tocase las latitudes de Hawái y el continente americano. 60 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos 4.2 Esquema del método y diferencias con el método existente El esquema general de nuestro método es parecido al de Muszynski, G. et al., primero realizamos análisis topológico de datos (ATD) y extraemos un vector por cada imagen con información topológica y luego clasificamos esos vectores con aprendizaje automático. Además, estamos realizando el método para el Pacífico noreste, pero este método es adaptable a otras regiones (al igual que con el método anterior). La diferencia reside en que empleamos otro método de ATD y de aprendizaje automático. El método consta de dos etapas: Etapa 1: Se extraen características topológicas de las imágenes con Análisis Topológico de Datos. Se emplea homología persistente en complejos cúbicos filtrados con el fin de extraer información sobre la evolución de las componentes conexas y los huecos. Estas características topológicas serán el input del clasificador de la segunda etapa. Etapa 2: Clasificamos las características previas con un método aprendizaje automático, los K Vecinos más Cercanos (KNearest Neighbours, KNN). La clasificación se divide en dos pasos: Primero se entrena el modelo para separar las características topológicas en dos clases (RA y no RA) en el conjunto de entrenamiento y después se le aplica el modelo al conjunto test. Posteriormente, se evalúa la clasificación en términos de exactitud, sensibilidad y precisión. Como estamos realizando aprendizaje supervisado, tanto para entrenar el modelo como para evaluar sus resultados, vamos a emplear las etiquetas que nos aporta el modelo TECA-BARD. El método tiene los siguientes datos de entrada y de salida. Datos de entrada: nimágenes en escalas de grises bidimensionales con datos en ninstantes de tiempo sobre la densidad superficial de vapor de agua (DSVA) en kgm−2. Cada imagen en escala de grises viene representada por una matriz. Para cada imagen tenemos una etiqueta (las obtenidas del modelo TECA-BARD) que nos indica a que clase pertenece (RA o no RA) esa observación. 4. método propuesto 61 Datos de salida. Etiquetas binarias. Si es RA,RA = 1, si no lo es, RA = 0. Así las tres diferencias de nuestro método con el de Muszynski, G. et al. son: 1. Se emplea homología persistente en dimensiones 0y1en vez de emplear el algoritmo Union-Find para rastrear la evolución de las componentes conexas en los subconjuntos de nivel superior. Con homología persistente en dimensiones 0y1se rastrea la evolución de las componentes conexas y de los huecos en los subconjuntos de nivel superior. El objetivo de este cambio es detectar, además de la evolución de las componentes conexas, posibles patrones de huecos en los subconjuntos de nivel superior de los RA que los diferencien. Además, con el algoritmo Union-Find se saca distinto tipo de información sobre las componentes conexas (tamaño de la componente conexa más vie- ja) que la información que se saca con homología persistente (información relacionada con el nacimiento y muerte de todas las componentes) 2. Se emplea KNN en vez de MVS como método de aprendizaje automático. Esto es debido a que KNN da resultados ligeramente mejores. 3. Se emplean etiquetas surgidas del modelo TECA-BARD y no del TECA. Esta diferencia simplemente es debida a las dificultades de acceso a los datos. Queda pendiente (al igual que en el método de Muszynski, G. et al.) probar el método con etiquetas provenientes de otras fuentes. 4.3 Primera etapa, análisis topológico de datos En esta etapa solo se emplean las imágenes sin las etiquetas (es decir, las coordenadas y el valor de DSVA). Vamos a intentar extraer información sobre la estructura topológica de los datos que sea relevante para que, en la siguiente etapa, el método de aprendizaje automático sepa diferenciar entre RA y no RA. Como hemos dicho anteriormente, vamos a usar homología persistente en complejos cúbicos filtrados para extraer información sobre la evolución de componentes conexas y de huecos en los datos. Una manera directa de hacer esto es modelar las imágenes en escala de grises con el valor de DSVA como complejos cúbicos filtrados (hemos explicado en la subsección 2.1.3 como se hace) y aplicar homología persistente en esos complejos cúbicos (explicada en la subsección 2.1.4). 62 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Sin embargo, no vamos a aplicar esto directamente sobre las imágenes con el valor de DSVA sino que vamos a modificar esas imágenes. Vamos a explicar primero lo que se va a aplicar sobre esas imágenes modificadas y luego podremos explicar (en la subsección 4.3.1) porque se ha aplicado sobre esas imágenes modificadas y no sobre las originales. Supongamos que tenemos esas imágenes ya modificadas respecto de las originales con el valor de DSVA. Como las imágenes modificadas no van a ser muy distintas, podemos pensar como que estas imágenes modificadas nos indican igualmente la concentración de humedad a lo largo de la región. Nos interesa saber como evoluciona la forma de las zonas de la región con una concentración de humedad mayor que cierto umbral a medida que va disminuyendo este umbral. Esto es debido a que los ríos atmosféricos, que son zonas de la superficie de la tierra con alta concentración de humedad, tienen una forma característica: son bandas alargadas y estrechas. Así, estamos rastreando qué forma tienen estas regiones (con concentración de humedad mayor que cierto umbral) y para qué umbral se obtiene cada forma (ya que no es lo mismo una banda fina de humedad con un umbral medio-bajo que con un umbral alto). La manera de medir como evoluciona la forma de estas zonas de la región es rastreando la evolución de las componentes conexas y huecos (de estas zonas). La manera de medir esta evolución es modelando estas imágenes modificadas como complejos cúbicos filtrados y aplicar homología persistente en ellos. La forma de modelar un complejo cúbico filtrado está explicada en la sección 2.1.3 y la homología persistente en complejos cúbicos en la sección 2.1.4. Con la homología persistente vamos a rastrear el nacimiento y muerte de clases de equivalencia en los grupos de homología de dimensión 1y0según vayamos moviendo el umbral en los subcomplejos de nivel superior. A partir de la homología persistente vamos a obtener dos diagramas de persistencia: El diagrama de persistencia de dimensión 0, con la información sobre el nacimiento y muerte de clases de equivalencia en los grupos de homología de dimensión 0. Equivalentemente, el diagrama de persistencia de dimensión 0nos da información sobre nacimiento y muerte de componentes conexas en los subcomplejos de nivel superior. El diagrama de persistencia de dimensión 1, con la información sobre el nacimiento y muerte de clases de equivalencia en los grupos de homología de 4. método propuesto 63 dimensión 1. Equivalentemente, el diagrama de persistencia de dimensión 0 nos da información sobre nacimiento y muerte de huecos en los subcomplejos de nivel superior. Estos diagramas de persistencia se han calculado en una implementación en ordenador, de manera que no va a haber puntos del diagrama con multiplicidad mayor que 1salvo que haya varias clases que nazcan en el mismo umbral. Observación 4.1.En la subsección 2.1.4 se han definido las filtraciones (def. 2.13) empleando subcomplejos de nivel inferior, mientras que nosotros queremos ver como evolucionan los subcomplejos de nivel superior (definiendo las filtraciones de esa manera). Existen dos soluciones posibles: La primera es que, para adaptarse a las definiciones que hemos hecho en el capítulo de conceptos previos, considerar, en vez de las matrices que representan las imágenes modificadas, las matrices multiplicadas por −1al modelar los complejos cúbicos filtrados. Así, al calcular los subcomplejos de nivel inferior respecto la función de filtrado asociada a la matriz multiplicada por −1estamos calculando los subcomplejos de nivel superior de respecto la función de filtrado asociada a la matriz original, ya que es claro que: f−1([r, +∞)) = (−f)−1((−∞, r]) . Esta es la solución que se emplea en la implementación a ordenador La segunda manera sería considerar las filtraciones en la definición 2.13 como subconjuntos de nivel superior. Además, para que estos subconjuntos de nivel superior sean subcomplejos deberíamos cambiar en la definición de función monótona (def. 2.9) f(P)≤f(Q)por f(P)≥f(Q)y en la construcción T (def. 2.12) el mínimo por el máximo al definir la función de filtrado en cubos de dimensión inferior a la máxima. Esta es la solución que se aplica cuando se hacen a mano ejemplos, como en el ejemplo de la subsección 4.3.2 (salvo el diagrama de persistencia). 4.3.1 Obtención de función de filtrado Veamos ahora por qué tenemos que aplicar homología persistente sobre esas matrices modificadas y no sobre las originales. 64 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos El motivo es porque la homología persistente no tiene en cuenta que los ríos atmosféricos tienen que ir desde las latitudes de Hawái hasta el continente americano. Estamos tratando de igual forma a una componente que va desde/hasta alguna de esas regiones que a una completamente alejada de ambas. Para cambiar esto, necesitamos no tener en cuenta en la información topológica extraída componentes conexas irrelevantes para la detección de RAs, como también se hizo en el método de Muszynski, G. et al. Para no tener en cuenta componentes conexas irrelevantes para la detección de RAs, se modifican las imágenes en escala de grises de manera que al aplicar homología persistente no entren componentes conexas que no toquen ninguna de las dos subregiones de interés (latitudes de Hawái o el continente americano). Como la conexión en grafos es equivalente a la conexión en complejos cúbicos y como existen muchos algoritmos implementados en ordenador para rastrear los elementos de las componentes conexas de un grafo, vamos a calcular las imágenes en escala de grises modificadas modelando la rejilla donde tenemos datos como un grafo de igual manera que se hizo en el método de Muszynski, G. et al.. Es decir, cada punto de la imagen (o elemento de la matriz) es un nodo del grafo que está unido por aristas a los puntos que están inmediatamente a su izquierda, derecha, arriba y abajo. Podemos ver una representación de grafo subyacente a la rejilla en la figura 3.2. En el algoritmo que presentamos en esta subsección, vamos a cambiar la matriz con los valores originales de DSVA por unos valores modificados donde no va a haber para ningún umbral componentes conexas en el subcomplejo cúbico de nivel superior (estricto) que no toquen los puntos de la rejilla que se corresponden con el continente americano o las latitudes de Hawái. Para ello, los valores modificados van a ser el primer umbral para el cual ese punto de la rejilla entra en el conjunto de nivel superior (estricto) en una componente conexa que toque algún punto correspondiente con América o las latitudes de Hawái. Concretamente, el algoritmo que aplicamos a cada imagen en escala de grises (es decir, una matriz) que tenemos en los datos 2es el siguiente. Sea Gel grafo subyacente a la rejilla. Sea Ala matriz con los valores de DSVA en los puntos de la rejilla N=n(i, j)|i∈ {1, . . . , n}, j ∈ {1, . . . , m}o, 2Dada una matriz A∈Mn,m(R), denotamos el elemento en la posición (i, j)como Aij . 4. método propuesto 65 siendo Nlos nodos de G. Sea AH ⊂Nlos puntos de la rejilla que se corresponden con el continente Americano o con las latitudes de Hawái. Fes la matriz que obtenemos como resultado, la imagen modificada. Sea p > 0un paso. 3 El input al algoritmo es A. Inicializamos la matriz F=Afunción de filtrado resultado y la matriz: indicador acumulado Indacum = 0n,m ∈Mn,m(R). Para cada umbral u∈U:= {60,60 −1p, 60 −2p, . . . , 0}4en orden descendente: 1. Inicializamos la matriz indicador Indi = 0n,m ∈Mn,m(R). 2. Calculamos las componentes conexas Comp del subgrafo Subgrf que tiene los nodos NSubgrf ={(i, j)∈N|Aij > u}y las aristas que tiene son las de Gque unen nodos de NSubgrf . 3. Determinamos Cs ={c∈Comp|AH ∩NSubgrf ∩c=∅} el conjunto de las componentes conexas de Subgrf que contengan algún elemento de AH ∩NSubgrf . 4. Para cada (i, j)∈N: a) Si (i, j)pertenece a un elemento de Cs entonces: Indiij = 1 e Indacumij =Indacumij + 1. b)•Si Indiij = 0 entonces Fij =u. •Si Indiij = 1 yIndacumij = 1 entonces Fij =Fij. •si Indiij = 1 yIndacumij = 1 entonces Fij =u. Proceso explicado Para cada punto (i, j)va a ocurrir el siguiente proceso a lo largo del bucle. En los primeros valores de U,Indiij =Indacumij = 0 porque ese nodo del grafo, o bien no está en el conjunto de nivel superior (estricto) NSubgrf o bien no pertenece a una componente conexa de Cs (las que tocan los nodos de AH). Para esos valores de Uvamos a asignarle ucomo valor de la función de filtrado (realmente es una cota superior porque en cada paso la vamos disminuyendo más). Para un único valor de U,Indiij =Indacumij = 1, este valor es el primer umbral en el que (i, j)∈A, i.e. (i, j)∈NSubgrf y(i, j)pertenece a una componente 3En nuestra implementación hemos considerado p= 0.25. 4Empezamos los umbrales en 60 porque el valor de DSVA está acotado por 60kgm−2(salvo valores atípicos) y los acabamos en 0porque DSVA es un variable no negativa. 72 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Figura 4.6: Esquema de la etapa 1del método. El esquema representa todo lo que se aplica para cada imagen en escala de grises del conjunto de datos. 4. método propuesto 73 En la imagen mostrada en la figura 4.7 no hay un RA, según las etiquetas provenientes del método TECA-BARD. En la figura 4.7, podemos ver la imagen en escala de grises del ejemplo y a su derecha podemos ver la imagen resultado de aplicarle a la imagen ejemplo el algoritmo para la obtención de una función de filtrado, La principal diferencia entre las dos imágenes es que en la modificada se han reducido los valores máximos, obteniendo un colorado más homogéneo en las zonas donde hay valores medio-altos. Finalmente, debajo de estas dos imágenes en la figura podemos ver los diagramas de persistencia de dimensión 0y1jun- tos. Estos diagramas son los que se obtienen al aplicarle a la imagen modificada homología persistente. En la figura 4.8, podemos ver a la izquierda los diagramas de persistencia transformados a coordenadas nacimiento-persistencia (arriba el de dimensión 0, arriba el de dimensión 1). A la derecha de los diagramas de persistencia transformados podemos ver sus imágenes de persistencia asociadas. En ellas se nota que los puntos con más persistencia (los que están más arriba) están mucho más ponderados en la imagen, ya que los píxeles que están más arriba necesitan que haya menos puntos en él (en el diagrama transformados) para tener valores altos. 4.4 Segunda etapa, método de los K Vecinos más Cercanos En esta segunda etapa, una vez obtenidos los vectores resultado de aplicar la etapa anterior a cada imagen en escala de grises del conjunto de datos, debemos clasificar estos vectores según las imágenes de la que provienen se correspondan o no a un río atmosférico. Para ello, vamos a emplear aprendizaje automático y por ello necesitamos usar las etiquetas aportadas por el modelo TECA-BARD (etiquetas binarias que nos indican si las imágenes se corresponden a un RA o no) para entrenar al modelo de aprendizaje automático que empleemos y para medir la calidad de la clasificación que haga el modelo. Por ello, el input a este método son: por un lado, los vectores resultado de aplicar la etapa anterior a cada imagen en escala de grises y por otro, las etiquetas binarias resultado de aplicar el modelo TECA-BARD a cada imagen del conjunto de datos. 74 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Figura 4.7: A la izquierda podemos ver una imagen en escala de grises donde se mide el valor de DSVA en kgm−2, esta imagen es una de las extraídas de nuestro conjunto de datos salidos del MERRA-2. A la derecha podemos ver la imagen en escala de grises resultado de aplicarle el algoritmo para la obtención de una función de filtrado a la imagen de la izquierda. Abajo, los diagramas de persistencia (dimensión 0en rojo, dimensión 1en azul) resultado de aplicarle a la imagen de la izquierda homología persistente en complejos cúbicos. 4. método propuesto 75 Figura 4.8: Arriba a la izquierda podemos ver el diagrama de persistencia de dimensión 0de la figura 4.7 transformado a coordenadas nacimiento-persistencia. Arriba a la derecha podemos ver la imagen de persistencia (para este ejemplo hemos considerado excepcionalmente un tamaño de píxel 1.75) del diagrama de dimensión 0de la figura 4.7. Abajo a la izquierda podemos ver el diagrama de persistencia de dimensión 1de la figura 4.7 transformado a coordenadas nacimientopersistencia. Abajo a la derecha podemos ver la imagen de persistencia (para este ejemplo hemos considerado excepcionalmente un tamaño de píxel 1.75) del diagrama de dimensión 1de la figura 4.7. 76 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Con esos datos se pone en marcha un método de clasificación supervisada donde las variables predictoras son las correspondientes a las características topológicas extraídas previamente y la variable a predecir (a clasificar, en este caso, al ser variable cualitativa) es la variable RA (valdría 1si en la observación hay un río atmosférico y −1en caso contrario). Para más información y detalle sobre los métodos de clasificación supervisada, ver sección 2.3. Se ha hecho una validación cruzada anidada LTO, procedimiento descrito en la subsección 2.3.2, para determinar cúal es el mejor modelo (entre un conjunto de modelos de aprendizaje automático candidatos) que clasifica mejor estos vectores provenientes de la primera etapa. El mejor modelo es el método de los K Vecinos Más Cercanos (KNearest Neighbours KNN). El método SVM con núcleo radial (considerado en el método de Muszynski, G.) clasificó ligeramente peor que el KNN (veremos los resultados en la sección 4.6) y no lo hemos considerado. El KNN predice la clase a la que pertenece un dato del conjunto test como la clase mayoritaria en los kdatos del conjunto de entrenamiento más cercanos (en el espacio dónde están las variables predictoras). Para más detalles sobre el KNN consultar la sección 2.10. El KNN depende de dos opciones. La primera es la distancia escogida para medir qué observaciones del conjunto de entrenamiento están más cerca. Hemos utilizado la distancia euclidiana. La segunda es escoger K, el número de vecinos considerados en el conjunto de entrenamiento. Para escoger este valor de K, entre un conjunto de valores candidatos, hemos empleado validación cruzada (no anidada). Este procedimiento y sus motivaciones están descritos en la sección 2.3.1. La importancia particular del valor de Kestá explicada en la sección 2.10. Una vez entrenado el KNN y habiendo escogido el Kóptimo, predecimos los valores para el conjunto test y calculamos las medidas de calidad de clasificación (descritas en 2.3.3). Observación 4.2.Al igual que otros métodos de aprendizaje automático, el KNN, tiene a sobreajustarse a la clase mayoritaria, por tanto, tiende a clasificar bastante bien los datos de la clase mayoritaria y mal los datos de la clase minoritaria. Para resolver este problema se ha aplicado remuestreo para obtener conjuntos de entrenamiento con clases equilibradas. El remuestreo se ha realizado antes de la división en conjunto de entrenamiento y test. Este remuestreo también se ha hecho en la validación cruzada anidada LTO previa en la que hemos seleccionado el mejor método de aprendizaje automático. 4. método propuesto 77 4.5 Implementación y código Hemos programado cada etapa de este método en dos lenguajes distintos. El código de implementación está disponible en: https://github.com/davafe/Metodo-deteccion-rios-atmosfericos. La primera etapa se ha programado en Python. Para leer los archivos netCDF con los datos obtenidos de MERRA-2 se ha empleado la librería netCDF4 ([26]). Hemos programado el algoritmo de obtención de función de filtrado (el algoritmo de la sección 4.3.1) y hemos utilizado el paquete cv2 ([27]) para calcular las componentes conexas de los subgrafos de nivel. En este algoritmo hemos considerado un paso entre umbrales de p= 0.25. Para modelar las imágenes en escala de grises como complejos cúbicos y computar la homología persistente de esos complejos hemos empleado el módulo Cubical Complex de gudhi ([28]). Finalmente, para calcular y representar las imágenes de persistencia hemos usado el módulo Persistence Imager de persim ([29]). La segunda etapa se ha programado en R. Destacar que hemos empleado la implementación de KNN del paquete class ([30]). El tiempo de computación de la primera etapa con paso p= 0.25 fue de alrededor de 30 horas. Con paso p= 1 fue de 8horas aproximadamente. Por lo tanto, el tiempo de computación depende mucho del paso escogido. En la segunda etapa, la computación duró alrededor de 1hora. Estos tiempos de computación se alcanzaron en un ordenador personal con procesador Intel Core i5. 4.6 Resultados Una vez aplicado el KNN a los vectores surgidos de aplanar y juntar las imágenes de persistencia, se mide la calidad de su clasificación en tres de las medidas de bondad que se han definido en la sección 2.3.3: exactitud, precisión y sensibilidad. Estas tres medidas se miden (valga la redundancia) en el conjunto test. En el cuadro 4.2 podemos ver la exactitud para los dos conjuntos de datos empleados (realmente son el mismo solo que en uno consideramos imágenes cada 3horas y en otro solo 1de esas 8imágenes diarias). Podemos ver como con datos 78 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Cuadro 4.2: Tabla con la exactitud en el conjunto test del KNN, para los distintos conjuntos de datos considerados. Mod. climático Resol. temporal Resol. espacial Exactitud MERRA-2 3h 50 km 0.8324 MERRA-2 24h 50 km 0.7209 Cuadro 4.3: Tabla con la precisión y sensibilidad en el conjunto test del KNN, para los distintos conjuntos de datos considerados. Mod. climático Resol. temporal Resol. espacial Precisión Sensibilidad MERRA-2 3h 50 km 0.7816 0.9026 MERRA-2 24h 50 km 0.7171 0.6728 cada 3horas obtenemos un resultado bastante bueno, con más del 83 % de imágenes correctamente. Sin embargo, si reducimos la precisión temporal, la exactitud se reduce bastante, más de un 10 %, obteniendo una exactitud mejorable. En el cuadro 4.3 podemos ver para los dos conjuntos de datos empleados la precisión y sensibilidad. Con precisión horaria de 3hobtenemos una magnífica sensibilidad, del 0.9026, es decir, que más del 90 % de los ríos atmosféricos se detectan. En este mismo conjunto de datos obtenemos una precisión no tan buena, del 0.7816, i.e., de las imágenes clasificadas como ríos atmosféricos, el 78.16 % son realmente ríos atmosféricos. Si reducimos la precisión horaria, la sensibilidad baja mucho, más de un 22 % y la precisión no baja tanto, algo más del 6,4 %. 4.6.1 Importancia de las dimensiones empleadas en homología y del método de aprendizaje automático escogido Una vez vistos estos resultados generales, cabe preguntarse dos cosas sobre la importancia de las novedades en nuestro método: ¿Si solo calculásemos la homología en dimensión 0, como en el método anterior, obtendríamos el mismo resultado? Dicho de otra manera, ¿incluir 4. método propuesto 79 Cuadro 4.4: Tabla con las distintas exactitudes, precisiones y sensibilidades obtenidas por dos distintos métodos (KNN y SVM con núcleo radial) y empleando las imágenes de persistencia de dimensiones 0 y/o 1. Método Dimensiones empleadas Exactitud Precisión Sensibilidad KNN 0,1 0.8324 0.7816 0.9022 KNN 0 0.8091 0.7643 0.8693 KNN 1 0.7717 0.7109 0.8811 SVM Radial 0,1 0.8310 0.8183 0.8320 SVM Radial 0 0.8162 0.7970 0.8264 SVM Radial 1 0.7183 0.6956 0.7334 la información topológica sobre homología en dimensión 1realmente mejora la clasificación? ¿Realmente utilizar KNN es mejor que emplear SVM con núcleo radial? Hemos computado nuestro método empleando en la clasificación solo imágenes de persistencia de dimensión 0y solo utilizando las imágenes de persistencia de dimensión 1. Además, hemos probado el método de aprendizaje automático SVM radial (explicado en la sección 2.4) sobre imágenes de dimensión 0y/o 1. La calidad de las clasificaciones de cada combinación de métodos y dimensiones empleadas se puede ver en el cuadro 4.4 en términos exactitud, precisión y sensibilidad. En primer lugar, fijándonos solamente en el método KNN, vemos como usar homología persistente en dimensión 0y1mejora en más de un 2,2 % la exactitud, en más de un 1,5 % la exactitud y en más de un 4 % la sensibilidad respecto a si solamente utilizásemos homología persistente en dimensión 0. Además, empleando homología persistente solo en dimensión 1obtenemos una buena exactitud, 0.7717, que si bien no mejora la también buena exactitud con solo dimensión 0, sí que obtenemos mejor sensibilidad con solamente dimensión 1que con solamente dimensión 0. Resumidamente, la homología persistente en dimensión 1aporta una mejora en la clasificación, especialmente en términos de sensibilidad (i.e., en clasificar correctamente las imágenes en escala de grises correspondientes a ríos atmosféricos) En segundo lugar, haciendo la comparativa KNN contra SVM núcleo radial, 80 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos podemos ver que KNN obtiene una exactitud ligeramente superior, empleando ambas dimensiones de homología, y una bastante superior, utilizando solo la dimensión 1de homología, que la que se obtiene con SVM radial. Sin embargo, si solo se usase la dimensión 0la SVM obtiene una exactitud algo superior. En términos de sensibilidad, el KNN es mejor (independientemente de que dimensiones se usen), sin embargo, la SVM es mejor en términos de exactitud (salvo si se emplea solamente dimensión 1). Concluimos que, de manera general, KNN obtiene mejores resultados que SVM, salvo que se use solo la dimensión 0. Finalmente, recordar que el modelo de aprendizaje automático está sesgado por las etiquetas empleadas provenientes de TECA-BARD, ya que hemos asumido que estas etiquetas son "verdad" y hemos entrenado y calificado nuestro modelo con ellas. Queda pendiente de comprobar y usar el método con distintos etiquetados provenientes de otros métodos. 5 Comparación de métodos: Estabilidad y resultados numéricos En este capítulo compararemos el método de Muszynski, G. et al. con el propuesto en este trabajo en términos de estabilidad y de resultados numéricos. En la sección 5.1 demostramos la estabilidad de nuestro método con un paso suficientemente pequeño y probamos la no estabilidad del método de Muszynski, G. et al. aunque el paso tienda a 0. En la sección 5.2 comparamos los resultados empíricos de ambos métodos. 5.1 Estabilidad La estabilidad es una propiedad importante que nos indica si el resultado de un algoritmo o método cambia mucho si el input cambia ligeramente. Que un método sea estable nos garantiza que una pequeña alteración en el input no va a provocar grandes cambios en el resultado. En esta sección compararemos la estabilidad que tienen la primera etapa, en la que se hace ATD, el método de Muszynski, G. et al. (descrito en el capítulo 3) y el método propuesto en este trabajo (descrito en el capítulo 4). Como ambos métodos tienen como input unas imágenes en escala de grises (que son matrices) y como salida un vector de Rk, recordamos las definiciones de las normas de matrices y vectores que vamos a usar en esta sección. 88 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos empleando la construcción T a las imágenes en escala de grises AyB, respectivamente. Se tiene que: W1,1(Dgmk(A), Dgmk(B)) ≤9||A−B||1. Este teorema es el que vamos a emplear en la demostración, ya que en nuestro método hemos usado subcomplejos de nivel superior. También vamos a usar el teorema siguiente sobre estabilidad de las imágenes de persistencia (explicadas en la sección 2.2): Teorema 5.4. Sean ByB′dos diagramas de persistencia y I(ρB)yI(ρB′) las imágenes de persistencia asociadas a ByB′, respectivamente, con función de distribución normal bivariante (con matriz de covarianzas diagonal y las varianzas de las dos variables valen σ2) y una función gde ponderación no negativa, continua, diferenciable a trozos y que vale cero en el eje horizontal del plano (si la persistencia es nula). Se tiene que: ||I(ρB)−I(ρB′)||∞≤5|∇g|+r10 π ||g||∞ σW1,∞(B, B′) donde ||g||∞= sup g(x, y)y|∇g|= sup ||∇g||2 Este teorema fue demostrado por Adams, H., Emerson, T en [15]. Con todos estos resultados podemos demostrar la estabilidad, cuando el paso tiende a 0, de la primera etapa de nuestro método: Teorema 5.5 (Estabilidad (con p→0) del método propuesto). Sean A, B ∈Mn,m(R)dos imágenes en escala de grises, p > 0un paso y h:C∈ Mn,m(R)7→ hC∈Rk, donde hCes el resultado de aplicarle a la matriz Cla primera etapa de nuestro método, la descrita en la sección 4.3. Se tiene que ||hA−hB||∞≤5 + r10 π609nm(||A−B||∞+p). Demostración. Sean A, B ∈Mn,m(R)yhA, hB∈Rklos resultados de aplicarle a las matrices AyB, respectivamente, el método descrito en la sección 4.3. Veamos paso por paso, de principio a final, todo el proceso que se les ha aplicado a las matrices y así vamos acotando en cada paso hasta acotar ||hA−hB||∞. 5. comparación de métodos: estabilidad y resultados numéricos 89 1. Lo primero que se les aplica a AyBes el algoritmo de obtención de función de filtrado (sec. 4.3.1) y, por el teorema 5.1: ||fA−fB||∞≤ ||A−B||∞+p(5.4) dónde fAyfBson las matrices resultantes de aplicar dicho algoritmo a A yB, respectivamente. 2. De fAyfBse obtienen Dgmk(A)yDgmk(B)con k= 0 yk= 1.Dgmk(A) yDgmk(B)son dos diagramas de persistencia de dimensión k∈ {0,1, . . . } asociados a aplicar homología persistente, con filtraciones de subcomplejos de nivel superior, a los complejos cúbicos filtrados asociados empleando la construcción T a las imágenes en escala de grises fAyfB. Podemos aplicar entonces el teorema 5.3 y obtener: W1,1(Dgm0(fA), Dgm0(fB)) ≤9||fA−fB||1. W1,1(Dgm1(fA), Dgm1(fB)) ≤9||fA−fB||1. Aplicando las desigualdades (5.2) y (5.4) obtenemos que: W1,1(Dgm0(fA), Dgm0(fB)) ≤9nm||fA−fB||∞≤9nm(||A−B||∞+p) (5.5) W1,1(Dgm1(fA), Dgm1(fB)) ≤9nm||fA−fB||∞≤9nm(||A−B||∞+p) (5.6) 3. De D0 A:= Dgm0(fA), D0 B:= Dgm0(fB)se obtienen I(ρD0 A)yI(ρD0 B).I(ρD0 A) yI(ρD0 B)son resultado de calcular las imágenes de persistencia a los diagramas de persistencia D0 A, D0 Ben las condiciones del teorema 5.4, con σ= 1 yg(x, y) = y, por tanto, por dicho teorema obtenemos que: ||I(ρD0 A)−I(ρD0 B)||∞≤5|∇g|+r10 π||g||∞W1,∞(D0 A, D0 B). De manera análoga, para D1 A:= Dgm0(fA), D1 B:= Dgm0(fB)se obtienen I(ρD1 A)yI(ρD1 B)cumpliendo que: ||I(ρD1 A)−I(ρD1 B)||∞≤5|∇g|+r10 π||g||∞W1,∞(D1 A, D1 B). Como cualquier matriz resultante del algoritmo de función de filtrado va a tomar valores en el conjunto de umbrales U={60,60−p, 60−2p, . . . , 0}(por el lema 5.1), todo diagrama de persistencia que obtengamos en el segundo paso va a obtener entre 0y60. Además, en las imágenes de persistencia (ver 90 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos la sección 2.2) solo se emplea la función de ponderación gsobre los diagramas en coordenadas nacimiento persistencia. Por tanto, podemos restringir la función ga puntos de R2con la segunda coordenada, la de persistencia, entre 0y60. Además, como en este trabajo consideramos g(x, y) = y(la proyección a la coordenada de persistencia), podemos acotar g, obteniendo ||g||∞≤60. Por otra parte, como g(x, y) = y,∇g(x, y) = (0,1) para todo (x, y), por lo tanto, |∇g|= 1. Con estas dos consideraciones obtenemos que: ||I(ρD0 A)−I(ρD0 B)||∞≤5 + r10 π60W1,∞(D0 A, D0 B), ||I(ρD1 A)−I(ρD1 B)||∞≤5 + r10 π60W1,∞(D1 A, D1 B). La desigualdad (5.3) implica que: ||I(ρD0 A)−I(ρD0 B)||∞≤5 + r10 π60W1,1(D0 A, D0 B), ||I(ρD1 A)−I(ρD1 B)||∞≤5 + r10 π60W1,1(D1 A, D1 B). Finalmente, empleando las desigualdades (5.5) y (5.6) obtenemos que: ||I(ρD0 A)−I(ρD0 B)||∞≤5 + r10 π609nm(||A−B||∞+p)(5.7) ||I(ρD1 A)−I(ρD1 B)||∞≤5 + r10 π609nm(||A−B||∞+p)(5.8) 4. Finalmente, se aplanan las matrices I(ρD0 A)yI(ρD1 A)y se empalman, los dos vectores resultantes obteniendo el vector resultado final hA. De manera análoga obtenemos el vector resultado final hBa partir de I(ρD0 B)yI(ρD1 B). Por las propiedades de la norma infinito de vectores y matrices es claro que: ||hA−hB||∞= m´ax{||I(ρD0 A)−I(ρD0 B)||∞,||I(ρD1 A)−I(ρD1 B)||∞}. De esta igualdad y de las desigualdades (5.7) y (5.8), concluimos que: ||hA−hB||∞≤5 + r10 π609nm(||A−B||∞+p) rematando la demostración del teorema. 5. comparación de métodos: estabilidad y resultados numéricos 91 Vamos a demostrar un teorema consecuencia del anterior, la estabilidad con paso fijo menor que la precisión de los datos. Teorema 5.6 (Estabilidad con paso menor que la precisión de los datos del método propuesto). Sean A, B ∈Mn,m(R)dos imágenes en escala de grises, p > 0un paso menor o igual que la precisión de los datos y h:C∈ Mn,m(R)7→ hC∈Rk, donde hCes el resultado de aplicarle a la matriz Cla primera etapa de nuestro método, la descrita en la sección 4.3. Se tiene que: ||hA−hB||∞≤5 + r10 π6018nm||A−B||∞. Demostración. Sean A, B ∈Mn,m(R)dos imágenes en escala de grises, p > 0un paso menor o igual que la precisión de los datos y h:C∈Mn,m(R)7→ hC∈Rk, donde hCes el resultado de aplicarle a la matriz Cla primera etapa de nuestro método. Si A=B,hA=hB, por tanto, la desigualdad se cumple trivialmente. Si A=B, como el paso p > 0es menor o igual que la precisión de los datos, entonces existe algún elemento de A−Bque es mayor que el paso p(ya que si no existiese, tendríamos que A=B). En consecuencia, ||A−B||∞≤p. Como cumplimos las condiciones del teorema 5.5 obtenemos que: ||hA−hB||∞≤5 + r10 π609nm(||A−B||∞+p) y como (||A−B||∞+p)≤2(||A−B||∞)obtenemos que se cumple la desigualdad a demostrar. 5.1.2 Carencia de estabilidad del método de Muszynski, G. et al. Para empezar, vamos a demostrar que la primera etapa del método de Muszynski, G. et al. no es estable cuando el paso tiende a cero. Es decir, vamos a demostrar que, por muy pequeño que sea el paso p, el algoritmo no va a ser estable. Para demostrarlo, es suficiente con ver que: ||A−B||∞→0, p →0⇏||FA−FB||∞→0(5.9) 92 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Figura 5.1: Grafo con el que modelamos las rejillas asociadas a las matrices. Las subregiones de interés están coloreadas (la zona naranja representa América y la zona violeta representa las latitudes de Hawái). donde FA, FB∈Rkson los vectores obtenidos de aplicarle esta primera etapa, con paso entre umbrales p > 0, a las matrices A, B ∈Mn,m(R), respectivamente. Por simplicidad, lo vamos a demostrar con matrices 5×5, pero es claramente extensible a las matrices de datos que estamos usando. Modelamos la rejilla asociada a las matrices como un grafo, considerando las subregiones de interés de la misma forma que hemos hecho en la subsección 3.3.2. Podemos ver este grafo en la figura 5.1. Consideremos una sucesión de valores δn=pn=1 n, con n∈N. Para demostrar (5.9), basta con encontrar para cada n∈Ndos matrices An, Bn∈M5,5(R)tal que ||An−Bn||∞< δnpero que ||FAn−FBn||∞≥Csiendo C > 0una constante. Remarcar que FAnyFBnson los vectores obtenidos de aplicarle esta primera etapa, con paso entre umbrales pn>0, a las matrices AnyBn, respectivamente. Consideremos, por tanto, las siguientes matrices: An=       10 10 10 10 10 10 10 10 10 10 10 10 25 + 1 4n25 25 10 25 25 10 10 25 25 10 10 10        , Bn=       10 10 10 10 10 10 10 10 10 10 10 10 25 −1 4n25 25 10 25 25 10 10 25 25 10 10 10        . 5. comparación de métodos: estabilidad y resultados numéricos 93 Figura 5.2: Grafos asociados a An(izquierda) y Bn(derecha). Las subregiones de interés están coloreadas (la zona naranja representa América y la zona violeta representa las latitudes de Hawái). Es claro que: An−Bn=       0 0 0 0 0 0 0 0 0 0 0 0 1 2n0 0 0 0 0 0 0 0 0 0 0 0        por lo que ||An−Bn||∞=1 2n< δn. Podemos ver los grafos asociados a AnyBnen la figura 5.2. Como el paso entre umbrales es p=1 n, el conjunto de umbrales considerados, U:= {60,60−p, 60−2p, . . . , 0}(en este orden), es tal que 25 ∈Uy que (25,25+ 1 4n]∩U=∅. Por tanto, resulta claro que el primer subgrafo de nivel superior1no vacío aparece, tanto para Ancomo Bn, con el umbral 25. Este subgrafo de nivel superior con umbral 25 lo podemos ver en la figura 5.3, tanto para Ancomo para Bn. Como 25 ∈U, va a existir una componente de los vectores FAnyFBnasociada al subgrafo de nivel superior de umbral 25, denotémosla por F25 AnyF25 Bn(para FAny 1Estamos considerando los subgrafos de nivel superior con ≥, no con >. 94 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos Figura 5.3: Subgrafos de nivel superior de umbral 25 asociados a An(izquierda) yBn(derecha). Las subregiones de interés están coloreadas (la zona naranja representa América y la zona violeta representa las latitudes de Hawái) FBn, respectivamente). Si demostramos que F25 An=7 25 yF25 Bn= 0 para todo n∈N, concluiríamos que ||FAn−FBn||∞≥C=7 25 >0y acabaríamos la demostración de inestabilidad. Veamos que F25 An=7 25. El primer umbral para el que el subgrafo de nivel superior tiene alguna componente conexa es 25, y esta componente toca ambas subregiones de interés (las que representan las latitudes de Hawái y el continente americano, coloreadas en figura 5.3), en consecuencia, si que la podemos tener en cuenta. Por lo tanto, como esta componente conexa es la que nace y podemos tener en cuenta antes, el valor de F25 Anes la proporción de nodos que tiene la componente conexa en el subgrafo de nivel (7) respeto del total de nodos del grafo (25). Veamos que F25 Bn= 0. Si bien existen dos componentes conexas en el subgrafo de nivel (ver figura 5.3), no podemos tener en cuenta ninguna de esas componentes conexas, ya que ninguna de ellas toca las dos subregiones de interés (las que representan las latitudes de Hawái y el continente americano, coloreadas en la figura). Por tanto, a los ojos del algoritmo, en el umbral 25 aún no apareció ninguna componente conexa a considerar, por lo que F25 Bn= 0. Observación 5.1.Es importante destacar que el método sería inestable aunque se tuviesen en cuenta las componentes conexas que no intersecan ambas subregiones de interés. Si para las mismas sucesiones δnypnconsiderásemos: 5. comparación de métodos: estabilidad y resultados numéricos 95 An=       10 10 10 10 10 10 10 10 10 10 10 10 25 + 1 4n25 + 1 4n25 + 1 4n 10 25 + 1 4n25 + 1 4n10 10 25 + 1 4n25 + 1 4n10 10 10        , Bn=       10 10 10 10 10 10 10 10 10 10 10 10 25 −1 4n25 −1 4n25 −1 4n 10 25 −1 4n25 −1 4n10 10 25 −1 4n25 −1 4n10 10 10        , tendríamos que An−Bn=       00000 00000 0 0 1 2n 1 2n 1 2n 01 2n 1 2n0 0 1 2n 1 2n0 0 0        .(5.10) Obtendríamos, de igual forma, que ||An−Bn||∞=1 2n< δn,F25 An=7 25 (ya que el subgrafo de nivel superior de umbral 25 sería el mismo) y que F25 Bn= 0 (ya que el subgrafo de nivel superior de umbral 25 sería vacío). Por lo cual, concluimos que también se puede demostrar la inestabilidad sin emplear la regla de que las componentes conexas que no toquen las subregiones de interés no cuentan. 5.2 Comparación de resultados numéricos Para comparar los resultados numéricos del método de Muszynski, G. et al. tenemos que ceñirnos al único conjunto de datos que hemos utilizado ambos: el obtenido del reanálisis MERRA-2 con precisión temporal de 3hy espacial de 50km. Es importante destacar que mientras los resultados Muszynski, G. et al. se obtuvieron empleando un etiquetado de las imágenes surgido de un método (TECA) distinto al nuestro (TECA-BARD). Estos etiquetados se usan para entrenar la máquina y medir la calidad en la etapa del aprendizaje automático de ambos 96 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos métodos, por tanto sesgan los modelos (por tanto, al ser distintos los etiquetados se sesgan de distinta forma). Igualmente, la calidad de las clasificaciones es comparable, aunque con ese matiz. Queda pendiente, por tanto, realizar la comparación con otros conjuntos de datos y habiendo usado con el mismo etiquetado. Una vez ya sabemos las condiciones de la comparación podemos confrontar resultados. Comparemos, pues, la calidad de la clasificación del modelo de Muszynski, G. et al. con la de nuestro método, con las distintas posibilidades que hemos considerado anteriormente (utilizando dimensiones de homología 0y/o 1y usando KNN o SVM con núcleo radial). En el cuadro 5.1 podemos ver la comparación en términos de exactitud en el conjunto test y en el cuadro 5.2 en términos de sensibilidad y precisión en el conjunto test. En ambos cuadros, como Muszynski, G. et al. rastrean con el Algoritmo Union-Find únicamente las componentes conexas (igual que la homología de dimensión 0), consideramos que se usa dimensión 0. El método propuesto en nuestro trabajo mejora en exactitud en hasta 4de las posibilidades (en todas las posibilidades en las que se emplea homología de dimensión 0) al método de Muszynski, G. et al.. La posibilidad en la que esa mejora es mejor es la posibilidad propuesta en el trabajo: utilizar homología persistente en ambas dimensiones y KNN. Con esta configuración, se mejora la exactitud en más del 3,2 %. El método propuesto en nuestro trabajo empeora en precisión en todas las posibilidades al método de Muszynski, G. et al., pero, sin embargo, lo mejora en sensibilidad en 5de las 6posibilidades. En particular, comparando solo con la posibilidad que hemos propuesto, nuestro método mejora en sensibilidad más de un 15 % y empeora en precisión más de un 5,8 % respeto del método de Muszynski, G. et al.. Resumidamente, el método propuesto (empleando KNN sobre imágenes de persistencia relativas a dimensiones 0y1) obtiene mejores resultados en exactitud (porcentaje de aciertos) y sensibilidad (porcentaje de aciertos en ríos atmosféricos). Sin embargo, obtiene peor resultado en precisión (porcentaje de ríos atmosféricos entre el total de imágenes clasificadas como río atmosférico). Por lo tanto, salvo que se tenga como objetivo principal mejorar la precisión, nuestro método obtiene mejor clasificación que el de Muszynski, G. et al.. 5. comparación de métodos: estabilidad y resultados numéricos 97 Cuadro 5.1: Tabla con las exactitudes (en el conjunto test) obtenidas por el método de Muszynski, G. et al. (el que utiliza en la primera etapa el Algoritmo Union-Find) y por nuestro método (empleamos homología persistente en la primera etapa) con las diversas posibilidades en las dimensiones consideradas en homología y el método en la segunda etapa. El método de Muszynski, G. et al. y el que proponemos nosotros (con KNN y usando ambas dimensiones en homología persistente) están destacados en negrita. Método 1ªetapa Método 2ªetapa Dim. usadas Exactitud Alg. Union-Find SVM Radial 0 0.80 Homología pers. KNN 0,1 0.8324 Homología pers. KNN 0 0.8091 Homología pers. KNN 1 0.7717 Homología pers. SVM Radial 0,1 0.8310 Homología pers. SVM Radial 0 0.8162 Homología pers. SVM Radial 1 0.7183 Cuadro 5.2: Tabla con las distintas precisiones y sensibilidades obtenidas (en el conjunto test) por el método de Muszynski, G. et al. (el que emplea en la primera etapa el Algoritmo Union-Find) y por nuestro método (empleamos homología persistente en la primera etapa) con las diversas posibilidades en las dimensiones consideradas en homología y el método en la segunda etapa. El método de Muszynski, G. et al. y el que proponemos nosotros (con KNN y empleando ambas dimensiones en homología persistente) están destacados en negrita. Método 1ªetapa Método 2ªetapa Dim. usadas Precisión Sensibilidad Alg. Union-Find SVM Radial 0 0.84 0.74 Hom. pers. KNN 0,1 0.7816 0.9022 Hom. pers. KNN 0 0.7643 0.8693 Hom. pers. KNN 1 0.7109 0.8811 Hom. pers. SVM radial 0,1 0.8183 0.8320 Hom. pers. SVM radial 0 0.7970 0.8264 Hom. pers. SVM radial 1 0.6956 0.7334 104 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos M., Sarangi, C., Sellars, S., Gorodetskaya, I., Kashinath, K., Kurlin, V., Mahoney, K., Muszynski, G., Pierce, R., Subramanian, A. C., Tome, R., Waliser, D., Walton, D., Wick, G., Wilson, A., Lavers, D., Prabhat, M., Collow, A., Krishnan, H., Magnusdottir, G., and Nguyen, P.: Atmospheric River Tracking Method Intercomparison Project (ARTMIP): project goals and experimental design, Geosci. Model Dev., 11, 2455-2474, https: //doi.org/10.5194/gmd-11-2455-2018, 2018 [9] Liu, Y., Racah, E., Prabhat, M., Correa, J., Khosrowshahi, A., Lavers, D., Kunkel, K., Wehner, M., and Collins, W. D.: Application of deep convolutional neural networks for detecting extreme weather in climate datasets, arXiv:1605.01156, 2016 [10] Tymochko, S., Munch, E., Dunion, J., Corbosiero, K. and Torn, R.: Using Persistent Homology to Quantify a Diurnal Cycle in Hurricane Felix Pattern Recognition Letters, 133, 137-143, 2020 [11] Choe, S. and Ramanna, S.: Cubical Homology-Based Machine Learning: An Application in Image Classification., Axioms, 11, 112. https://doi.org/10.3390/axioms11030112, 2022 [12] Kaczynski, T., Mischaikow, K. and Mrozek, M.: Computational Homology , volume 157 of Aplied Mathematical Sciences. Springer-Verlag, 2004 [13] Bleile, B., Garin, A., Heiss, T., Maggs, K. and Robins, V.: The Persistent Homology of Dual Digital Image Constructions., arXiv:2102.11397, https: //doi.org/10.48550/arXiv.2102.11397, 2022 [14] Edelsbrunner H. and Harer, J.: Computational Topology, An Introduction. American Mathematical Society, 2010. [15] Adams, H., Emerson, T., Kirby, M., Neville, R., Peterson, C., Shipman, P., Chepushtanova, S., Hanson, E., Motta, F. and Ziegelmeier, L.: Persistence Images: A Stable Vector Representation of Persistent Homology, Journal of Machine Learning Research, 18, 16-337, http://jmlr.org/papers/v18/ 16-337.html, 2017 [16] Russell, S. and Norvig, P.: Artificial Intelligence: A Modern Approach. Pearson, London, 2021. BIBLIOGRAFÍA 105 [17] Dinh, T. L. A. and Aires, F.: Nested leave-two-out cross-validation for the optimal crop yield model selection, Geosci. Model Dev., 15, 3519–3535, https://doi.org/10.5194/gmd-15-3519-2022, 2022 [18] James, G., Witten, D., Hastie, T., Tibshirani, R.: An Introduction to Statistical Learning : with Applications in R. Springer, New York, 2013. [19] Hastie, T., Tibshirani, R., and Friedman, J.: The elements of statistical learning: Data Mining, Inference and Prediction. Springer Verlag, Berlín, 2009 [20] Eaton, B.: User’s Guide to the Community Atmosphere Model CAM-5.1, NCAR, 2011. [21] Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.- K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, S., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The modern-era retrospective analysis for research and applications, version 2 (MERRA-2), J. Climate, 30, 5419- 5454, 2017 [22] Prabhat, M., Byna, S., Vishwanath, V., Dart, E., Wehner, M., and Collins, W. D.: TECA: Petascale pattern recognition for climate science, in: International Conference on Computer Analysis of Images and Patterns, Springer, 426-436, 2015 [23] Muszynski, G., Kashinath, K., Kurlin, V., Wehner, M., and Prabhat, M.: NERSC Science Gateway,https://doi.org/10.25342/GMD_2018, 2018 [24] Hopcroft, J. E. and Ullman, J. D.: Set merging algorithms, SIAM J. Comput., 2, 294-303, 1973 [25] O’Brien, T. A., Risser, M. D., Loring, B., Elbashandy, A. A., Krishnan, H., Johnson, J., Patricola, C. M., O’Brien, J. P., Mahesh, A., Prabhat, M., Arriaga Ramirez, S., Rhoades, A. M., Charn, A., Inda Díaz, H., and Collins, W. D.: Detection of atmospheric rivers with inline uncertainty quantification: TECA-BARD v1.0.1, Geosci. Model Dev., 13, 6131–6148, https://doi.org/10.5194/gmd-13-6131-2020, 2020 106 análisis topológico de datos y aprendizaje automático en la detección de ríos atmosféricos [26] Unidata, Network Common Data Form (NetCDF) 1.5.7 [software]. UCAR/Unidata Program Center, Boulder, EE.UU., https://doi.org/10. 5065/D6H70CW6, 2021 [27] Bradski, G.: The OpenCV Library. Journal of Software Tools, 120, 122-125, 2000 [28] Dlotko, P.: Cubical complex, GUDHI User and Reference, Manual, GUDHI Editorial Board„ http://gudhi.gforge.inria.fr/doc/latest/ group__cubical__complex.html, 2015 [29] Saul, N. and Tralie, C.: Scikit-TDA: Analysis for Python, Zenodo.: http: //doi.org/10.5281/zenodo.2533369, 2019 [30] Venables, W.N. and Ripley, B.D.: Modern Applied Statistics with S, Cuarta Edición, Springer, New York, ISBN 0-387-95457-0 https://www.stats.ox. ac.uk/pub/MASS4/, 2002 [31] Skraba, P. and Turner, K.: Wasserstein Stability for Persistence Diagrams, arXiv: Algebraic Topology, https://doi.org/10.48550/arXiv. 2006.16824, 2020