Câmpus de São José do Rio Preto Bianca Barreiro Clusterings Hierárquicos em Networks e Aplicações São José do Rio Preto 2019 Bianca Barreiro Clusterings Hierárquicos em Networks e Aplicações Dissertação apresentada como parte dos requisitos para obtenção do título de Mestre em Matemática, junto ao Programa de Pós- Graduação em Matemática, do Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Orientador: Prof. Dr. Thiago de Melo Financiadora: CAPES São José do Rio Preto 2019 B271c Barreiro, Bianca Clusterings Hierárquicos em Networks e Aplicações / Bianca Barreiro. -- São José do Rio Preto, 2019 82 p. : il., tabs. Dissertação (mestrado) - Universidade Estadual Paulista (Unesp), Instituto de Biociências Letras e Ciências Exatas, São José do Rio Preto Orientador: Thiago de Melo 1. Análise Topológica de Dados. 2. Clustering Hierárquico. 3. Networks. I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca do Instituto de Biociências Letras e Ciências Exatas, São José do Rio Preto. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. Bianca Barreiro Clusterings Hierárquicos em Networks e Aplicações Dissertação apresentada como parte dos requisitos para obtenção do título de Mestre em Matemática, junto ao Programa de Pós- Graduação em Matemática, do Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Financiadora: CAPES Comissão Examinadora Prof. Dr. Thiago de Melo Orientador Sergio Tsuyoshi Ura UNESP - Câmpus de Rio Claro Marcio Fuzeto Gameiro USP - Câmpus de São Carlos Rio Claro 18 de fevereiro de 2019 AGRADECIMENTOS Quero dedicar esse espaço para agradecer a todas as pessoas que foram im- portantes para mim desde o início da minha jornada até o presente momento. Quero agradecer aos meus pais, Luiz Antonio Barreiro e Lucimara Cristina Fo- saluza Barreiro, que sempre me apoiaram e acreditaram em minha capacidade e que fizeram inúmeros sacrifícios para que eu pudesse chegar onde estou. Ao meu irmão Felipe, que mesmo ainda pequeno, me fez amadurecer e criar responsabilidades. Ao meu namorado e melhor amigo Levi, que sempre me incentivou a continuar nessa jornada, nunca duvidou do meu potencial, até mesmo quando eu duvidei. Ao meu orientador, Prof. Dr. Thiago de Melo, que teve paciência pra me ensinar sobre programação, acreditou em meu potencial e confiou em meu trabalho. À Profa. Dra. Alice Kimie Miwa Libardi, que foi minha primeira orientadora, me fez apreciar a beleza de matemática e acreditou na minha capacidade. A todos os outros professores que passaram pela minha graduação e pós- graduação e que de alguma forma fizeram com que eu me apaixonasse cada vez mais pela matemática e pelo compromisso de sempre querer aprender coisas no- vas. Dedico um agradecimento especial a todos os meus amigos de graduação e pós-graduação, que me ajudaram e me acompanharam em toda esta jornada. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoa- mento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001. RESUMO Neste trabalho estudamos networks e métodos de Clustering Hierárquico de forma axiomática. Apresentamos alguns programas na linguagem Python aplicados à análise de dados de migração populacional, com o intuito de ilustrar os métodos estudados. Palavras-chave: Análise Topológica de Dados, Clustering Hierárquico, Networks. ABSTRACT In this work we study networks and an axiomatic construction of Hierarchical Clus- tering. We present some Python programs used to analyze human migration data and illustrate the studied methods. Keywords: Topological Data Analysis, Hierarchical Clustering, Networks. Lista de Figuras 1.1 A menor network não trivial . . . . . . . . . . . . . . . . . . . . . . . . 21 1.2 Network assimétrica de três nós . . . . . . . . . . . . . . . . . . . . . . 22 1.3 Network assimétrica de quatro nós . . . . . . . . . . . . . . . . . . . . 24 1.4 Exemplo de dendrogramas . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.1 Network canônica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Representação de um clustering recíproco . . . . . . . . . . . . . . . . . 39 3.2 Representação de um clustering não-recíproco . . . . . . . . . . . . . . 40 5.1 Representação da função top_migration_by_state() . . . . . . . . . 66 A.1 Fluxo migratório do estado de Nova Iorque . . . . . . . . . . . . . . . . 71 A.2 Dendrograma (parcial) de migração nacional dos EUA, método recí- proco, função peso f2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.3 Migração mundial, método recíproco, função peso f1 . . . . . . . . . . 73 A.4 Migração mundial, método não-recíproco, função peso f1 . . . . . . . . 74 A.5 Migração mundial, método recíproco, função peso f2 . . . . . . . . . . 75 A.6 Migração mundial, método não-recíproco, função peso f2 . . . . . . . . 76 Lista de Tabelas 5.1 Conteúdo de prepare_data_states_num_counties.csv . . . . . . . . 48 B.1 Lista de países, códigos ISO e continentes . . . . . . . . . . . . . . . . . 79 C.1 Lista de dez países com maior �uxo migratório . . . . . . . . . . . . . . 81 C.2 Códigos ISO e cores das regiões continentais . . . . . . . . . . . . . . . 81 C.3 Lista de dez condados com maior �uxo migratório . . . . . . . . . . . . 82 Sumário Introdução 19 1 Preliminares 21 2 Axiomas do Valor e da Transformação 27 2.1 Algumas equivalências . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3 Clustering Recíproco e Não-Recíproco 37 4 Algoritmos 43 4.1 Álgebra Dioide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Algoritmo para ultramétricas . . . . . . . . . . . . . . . . . . . . . . . 44 5 Aplicações 47 5.1 Migração da população dos EUA . . . . . . . . . . . . . . . . . . . . . 47 5.1.1 Preparando os dados: prepare_data.py . . . . . . . . . . . . . 49 5.1.2 Criando as ultramétricas: compute_ultrametric.py . . . . . . 52 5.1.3 Criando os dendrogramas: plot_dendrogram.py . . . . . . . . . 55 5.2 Migração da população mundial . . . . . . . . . . . . . . . . . . . . . . 58 5.2.1 Resumo dos dados iniciais . . . . . . . . . . . . . . . . . . . . . 58 5.2.2 Programa world-migration.py . . . . . . . . . . . . . . . . . . 58 5.3 Ranking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.3.1 Programa ranking_world_migration.py . . . . . . . . . . . . 62 5.3.2 Classi�cando os dados: ranking_USA_migration.py . . . . . . 65 Referências 69 A Dendrogramas 71 B Países, códigos ISO e continentes 77 C Ranking de migração 81 C.1 Migração mundial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 C.2 Migração dos EUA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Introdução Agrupar, armazenar e analisar dados é de fundamental importância em todas as áreas da Ciência. Quase todo tipo de informação pode ser transformada em um con- junto de dados. Por exemplo, a ideia de usar um conjunto de pontos (ou nuvem de pontos) para visualizar um objeto foi introduzida recentemente, mas não tornou-se popular até o surgimento de uma grande quantia de dados disponíveis para estudo. Portanto, o desenvolvimento de métodos para entendermos dados tão complexos se faz necessário, e um grande esforço para adaptar métodos topológicos a várias aplicações tem sido feito recentemente. Esta nova linha de pesquisa é chamada de Análise Topológica de Dados�TDA (do inglês Topological Data Analysis). Um conjunto de dados admite diferentes interpretações que dependem se esse con- junto de dados é métrico; simétrico, mas não necessariamente métrico, ou até mesmo assimétrico. Esses conjuntos de dados são denominados networks. Teoricamente, uma network é um conjunto �nito de nós juntamente com uma função dissimilaridade, em que cada nó está ligado a outro e existe um peso entre eles. Um nó está sempre ligado a ele mesmo, porém esse peso é nulo. Além disso, quando falamos de dois nós distintos x e x′, seus pesos não são necessariamente os mesmos, ou seja, o peso de x estar ligado a x′ pode ser diferente do peso de x′ a x. Quando temos um conjunto de dados, precisamos saber como interpretá-lo. Um dos meios para fazê-lo é construir métodos de Clustering Hierárquicos. Esses métodos servem para nos dizer quando os nós da network se agrupam mais rapidamente que outros, ou seja, quando temos um conjunto de dados de um problema real, dois dados se agruparem mais rapidamente signi�ca que eles têm mais a�nidade. Para interpretar de forma mais clara o que esses métodos nos dizem, usamos dendrogramas. Os dendrogramas podem ser interpretados como uma estrutura que produz dife- rentes clusterings em diferentes parâmetros. No Capítulo 5, usamos dados reais de migração dos EUA e migração mundial para produzir os dendrogramas descritos no Apêndice A. Neste trabalho, estudamos networks e métodos de Clustering Hierárquico de forma axiomática. Apresentamos alguns programas na linguagem Python aplicados à análise de dados de migração populacional, com o intuito de ilustrar os métodos estudados. A principal referência para esse projeto é [1]. 19 1 Preliminares Neste capítulo, apresentamos os principais conceitos sobre networks, necessários para as aplicações. Seja X um conjunto �nito não-vazio. Uma função AX : X×X → R+ é uma função dissimilaridade se AX(x, x′) = 0⇔ x = x′, ∀x, x′ ∈ X. Observação 1.1. Uma função dissimilaridade AX não é, necessariamente, uma métrica no conjunto �nito X, pois não se exige a validade da desigualdade triangular. Mais ainda, a função pode ser assimétrica, ou seja, AX(x, x′) 6= AX(x′, x), para algum x, x′ ∈ X. De�nição 1.2. Uma network NX é um par (X,AX) onde X é um conjunto �nito não-vazio e AX : X ×X → R+ é uma função dissimilaridade. O conjunto de todas as networks será denotado por N . Em alguns casos é conveniente escrever X = {x1, x2, . . . , xn} e chamar cada ele- mento de nó, além de interpretar a função dissimilaridade AX como uma matriz AX ∈ M(n × n,R+), com (AX)i,j = AX(xi, xj), para 1 ≤ i, j ≤ n. Note que os elementos da diagonal (AX)i,i = AX(xi, xi) são zero. Networks em N podem ter diferentes conjuntos de nós X, assim como diferentes funções dissimilaridades AX . Se X = {x1}, necessariamente AX é nula e, neste caso, a network (X,AX) é cha- mada trivial. As menores networks não triviais contêm dois nós p e q e duas dissimi- laridades α e β, como na Figura 1.1. p q α β Figura 1.1: A menor network não trivial De�nição 1.3. Considere a função dissimilaridade Ap,q dada por Ap,q(p, q) = α e Ap,q(q, p) = β, para algum α, β > 0 e de�na a network de dois nós −→ ∆2(α, β) com parâmetros α e β por −→ ∆2(α, β) := ({p, q}, Ap,q). (1.1) 21 22 Preliminares Exemplo 1.4. Considere a network da Figura 1.2, cujo conjunto de nós éX = {a, b, c}. Assim, podemos interpretar a função dissimilaridade como uma matriz (não-simétrica) AX = 0 2 5 4 0 1 3 7 0  . a b c 2 54 1 3 7 Figura 1.2: Network assimétrica de três nós De�nição 1.5. Um clustering do conjunto X é uma partição PX de X, isto é, uma coleção PX = {B1, . . . , BJ} de subconjuntos de X que são dois a dois disjuntos, isto é, Bi ∩ Bj = ∅ para i 6= j, e cobrem X, ou seja, ∪Ji=1Bi = X. Os conjuntos B1, . . . , BJ são chamados de clusters de PX . Uma partição PX = {B1, . . . , BJ} de X sempre induz, e é induzida por, uma relação de equivalência ∼PX em X, onde para todo x, x′ ∈ X temos que x ∼PX x′ se, e somente se, x e x′ pertencem a um mesmo cluster Bi, para algum i. Neste trabalho, vamos focar em métodos de clustering hierárquicos. O resultado desse tipo de método não é uma única partição PX , mas sim uma coleção DX de partições DX(δ) de X indexadas por um parâmetro δ ≥ 0, para a qual escrevemos x ∼DX(δ) x ′ se, e somente se, x e x′ estão em um mesmo cluster de DX(δ). Uma coleção DX é chamada de dendrograma se possui as seguintes propriedades: (D1) Condições de fronteira: Para δ = 0, a partição DX(0) possui apenas nós unitários de X, ou seja, DX(0) = {{x}, x ∈ X} e para algum δ0 su�cientemente grande, DX(δ0) possui todos os elementos de X em um único conjunto, isto é, DX(δ0) = {X}. (D2) Hierarquia: À medida que δ aumenta, os clusters podem ser combinados, mas nunca separados. Mais precisamente, para quaisquer δ1 < δ2, se x ∼DX(δ1) x ′ então x ∼DX(δ2) x ′, para quaisquer pontos x, x′ ∈ X. (D3) Continuidade à direita: Para todo δ ≥ 0, existe ε > 0 tal que DX(δ) = DX(δ′), ∀δ′ ∈ [δ, δ + ε]. Preliminares 23 Observação 1.6. Note que se o processo de clustering for via bola fechada, ou seja, x ∼dX(δ) x ′ ⇔ dX(x, x′) ≤ δ, em que dX é uma métrica, a condição (D1) é trivial. A segunda parte da condição (D1) juntamente com (D2) implica que devemos ter DX(δ) = {X} para todo δ ≥ δ0. Denotamos por [x]δ a classe de equivalência de x ∈ X no parâmetro δ, ou seja, [x]δ := {x′ ∈ X | x ∼DX(δ) x ′}. Por (D1), temos que [x]0 = {x} e [x]δ0 = X, para cada x ∈ X. Denotando por D o conjunto de todos os dendrogramas, um método de Clustering Hierárquico é uma função H : N → D, tal que o conjunto X é preservado, ou seja, a imagem de uma network sobre X é um dendrograma sobre X. Dada uma network NX = (X,AX), denotamos por DX = H(X,AX) sua imagem pelo método H. De�nição 1.7. Sejam (X,AX) uma network e x, x′ ∈ X. Uma cadeia de x a x′ é uma sequência ordenada de nós em X, escrita como C(x, x′) := [x = x0, x1, . . . , xl = x′], a qual dizemos que conecta x a x′. Os links da cadeia são as arestas que conectam nós consecutivos na direção imposta pela cadeia. De�nição 1.8. Dadas duas cadeias C(x, x′) = [x = x0, x1, . . . , xl = x′] e C(x′, x′′) = [x′ = x′0, x ′ 1, . . . , x ′ l′ = x′′], de�nimos a cadeia concatenada, C(x, x′) ] C(x′, x′′) := [x = x0, . . . , xl = x′ = x′0, . . . , x ′ l′ = x′′]. A operação ] de concatenação de cadeias é associativa, isto é, [C(x, x′) ] C(x′, x′′)] ] C(x′′, x′′′) = C(x, x′) ] [C(x′, x′′) ] C(x′′, x′′′)]. Além disso, observe que as cadeias C(x, x′) = [x = x0, x1, . . . , xl = x′] e sua inversa C(x′, x) := [x′ = xl, xl−1, . . . , x0 = x] são diferentes. De�nição 1.9. Dada uma cadeia C(x, x′) = [x = x0, x1, . . . , xl = x′], de�nimos o custo da cadeia como sendo max i|xi∈C(x,x′) AX(xi, xi+1). O custo mínimo de cadeias entre x e x′ é de�nido por ũ∗X(x, x′) := min C(x,x′) max i|xi∈C(x,x′) AX(xi, xi+1), (1.2) ou seja, é o custo mínimo entre todas as cadeias que conectam x a x′. O uso da notação ũ∗X se deve ao fato de que no Capítulo 2 veremos outro conceito que utiliza uX como notação. Em networks assimétricas, os custos mínimos de cadeia ũ∗X(x, x′) e ũ∗X(x′, x) em geral são diferentes, mas são iguais em networks simétricas. Nesse caso, os custos 24 Preliminares ũ∗X(x, x′) = ũ∗X(x′, x) fazem parte da de�nição do dendrograma single-linkage, que veremos mais adiante. De�nimos um ciclo como sendo uma cadeia da forma C(x, x), para algum x ∈ X, tal que C(x, x) contenha pelo menos um nó além de x. Como um ciclo é uma cadeia particular, o custo do ciclo é dado na De�nição 1.9. Além disso, de�nimos o custo mínimo de ciclos da network (X,AX) por mlc(X,AX) := min x min C(x,x) max i|xi∈C(x,x) AX(xi, xi+1). (1.3) De�nição 1.10. A separação de uma network (X,AX) é de�nida por sep(X,AX) := min x6=x′ AX(x, x′), (1.4) ou seja, é a menor dissimilaridade positiva. Note que, por (1.2) e (1.4), temos ũ∗X(x, x′) ≥ sep(X,AX). Assim, por (1.3), sep(X,AX) ≤ mlc(X,AX). (1.5) Além disso, no caso de networks simétricas, temos a igualdade em (1.5), pois mlc(X,AX) = min x min C(x,x) max i|xi∈C(x,x) AX(xi, xi+1) = min x min C(x,x) max i|xi∈C(x,x) AX(xi+1, xi) ≤ min C(x,x′)]C(x′,x) x6=x′ maxAX(xi, xi+1) ≤ min x 6=x′ AX(x, x′) = sep(X,AX). Exemplo 1.11. Considere a network NX = (X,AX), com conjunto de nós X = {a, b, c, d}, representada pela Figura 1.3. a b c d 2 5 3 1 4 9 7 10 4 3 2 8 Figura 1.3: Network assimétrica de quatro nós A função dissimilaridade é representada pela matriz AX =  0 2 5 3 4 0 9 1 10 4 0 7 8 2 3 0  , Preliminares 25 da qual vemos que sep(X,AX) = 1. As Figuras 1.4a e 1.4b representam os dendrogramas recíproco e não-recíproco (que veremos mais à frente) que podem ser derivados da network acima. a bc d 0 1 2 3 4 5 6 7 8 (a) Dendrograma recíproco a bc d 0 .5 1 1.5 2 2.5 3 3.5 4 (b) Dendrograma não-recíproco Figura 1.4: Exemplo de dendrogramas 2 Axiomas do Valor e da Transformação Neste capítulo, estudaremos os Axiomas do Valor e da Transformação e ultramé- tricas, que se fazem necessários para podermos estudar métodos de Clustering Hierár- quico. Lembremos da de�nição de −→ ∆2(α, β) vista em (1.1), −→ ∆2(α, β) = ({p, q}, Ap,q), em que Ap,q(p, q) = α e Ap,q(q, p) = β. De�nição 2.1. Dadas duas networks NX = (X,AX) e NY = (Y,AY ), uma função φ : X → Y é chamada redução de dissimilaridade se AX(x, x′) ≥ AY (φ(x), φ(x′)), ∀x, x′ ∈ X. (A1) Axioma do Valor. Dizemos que H satisfaz o Axioma do Valor se o dendrograma Dp,q = H( −→ ∆2(α, β)) produzido por H aplicado à network −→ ∆2(α, β) satisfaz Dp,q(δ) = { {{p}, {q}}, para 0 ≤ δ < max(α, β), {{p, q}}, para δ ≥ max(α, β). (A2) Axioma da Transformação. Considere duas networks NX = (X,AX) e NY = (Y,AY ) e uma função redução de dissimilaridade φ : X → Y . Dizemos que o método H satisfaz o Axioma da Transformação se DX = H((X,AX)) e DY = H((Y,AY )) são dendrogramas tais que se x ∼DX(δ) x ′ então φ(x) ∼DY (δ) φ(x′), com δ ≥ 0. Dizemos que um método de clustering hierárquico H é admissível se valem os axiomas (A1) e (A2) para H. De�nição 2.2. Seja X um conjunto de nós. Uma ultramétrica em X é uma função não-negativa uX : X ×X → R+ satisfazendo: (i) Identidade. uX(x, x′) = 0 se, e somente se, x = x′, para todo x, x′ ∈ X. (ii) Simetria. uX(x, x′) = uX(x′, x), para todo x, x′ ∈ X. (iii) Desigualdade triangular forte. Para quaisquer x, x′, x′′ ∈ X, uX(x, x′′) ≤ max(uX(x, x′), uX(x′, x′′)). (2.1) 27 28 Axiomas do Valor e da Transformação Como (2.1) implica que uX(x, x′′) ≤ uX(x, x′)+uX(x′, x′′), para quaisquer x, x′, x′′ ∈ X, segue que conjuntos ultramétricos são casos particulares de conjuntos métricos. O interesse principal em estudar ultramétricas é estabelecer uma estrutura que preserva bijeção entre o conjunto dos dendrogramas e o conjunto das ultramétricas. A seguir, de�niremos uma função Ψ : D → U do conjunto D dos dendrogramas no conjunto U das networks dotadas com ultramétricas. Mais precisamente, dado um dendrograma DX , seja Ψ(DX) = (X, uX), onde a ultramétrica uX é dada por uX(x, x′) := min{δ ≥ 0 | x ∼DX(δ) x ′}. (2.2) Por outro lado, de�niremos também uma função Υ : U → D. Dada uma ultra- métrica uX no conjunto �nito X e dado δ ≥ 0, seja a relação ∼uX(δ) em X dada por x ∼uX(δ) x ′ ⇔ uX(x, x′) ≤ δ. (2.3) Daí, tomamos DX(δ) := Xmod ∼uX(δ) e de�nimos Υ(X, uX) := DX . Isto posto, tem-se: Teorema 2.3. As funções Ψ : D → U e Υ : U → D estão bem de�nidas. Além disso, Ψ ◦Υ é a identidade em U e Υ ◦Ψ é a identidade em D. Demonstração. Será feita em duas partes. Parte 1. Primeiramente, mostremos que Ψ está bem de�nida. Para isso, basta mostrar que uX é uma ultramétrica. Notemos que o mínimo em (2.2) existe. De fato, seja γ = inf{δ ≥ 0 | x ∼DX(δ) x ′} e suponhamos que γ /∈ {δ ≥ 0 | x ∼DX(δ) x ′}. Por (D3) aplicado a γ temos que existe ε > 0 tal que DX(δ) = DX(γ), para todo δ ∈ [γ, γ + ε]. Em particular, DX(γ) = DX(γ + ε), implicando que γ + ε é cota inferior do conjunto acima, o que contradiz o fato de γ ser o ín�mo de tal conjunto. Portanto, γ é elemento mínimo. Agora, provemos as condições (i)-(iii) da De�nição 2.2. (i) Temos que uX(x, x′) = 0⇔ x ∼DX(0) x ′ ⇔ x = x′, pois ∼DX(δ) é uma relação de equivalência. (ii) Como ∼DX(δ) é uma relação de equivalência, temos que x ∼DX(δ) x ′ ⇔ x′ ∼DX(δ) x. Logo, uX(x, x′) = uX(x′, x). (iii) Dados x, x′, x′′ ∈ X, consideremos uX(x, x′) = δ1 e uX(x′, x′′) = δ2. Tomando δ = max{δ1, δ2}, temos que x ∼DX(δ1) x ′ implica x ∼DX(δ) x ′ e que x′ ∼DX(δ2) x ′′ implica x′ ∼DX(δ) x ′′. Pela transitividade da relação ∼DX(δ), segue que x ∼DX(δ) x′′, ou seja, uX(x, x′′) ≤ δ. Logo, uX é uma ultramétrica em X e Ψ é bem de�nida. A seguir, para provarmos que Υ está bem de�nida, primeiro mostremos que ∼uX(δ) dada em (2.3) é uma relação de equivalência, e portanto induz uma partição do conjunto X, para cada δ ≥ 0. � ∼uX(δ) é re�exiva, ou seja, x ∼uX(δ) x, pois 0 = uX(x, x) ≤ δ, ∀δ ≥ 0. � ∼uX(δ) é simétrica, pois dados x, x′ ∈ X, se x ∼uX(δ) x ′ então uX(x, x′) ≤ δ, e como uX é ultramétrica, uX(x, x′) = uX(x′, x). Portanto x′ ∼uX(δ) x. Axiomas do Valor e da Transformação 29 � ∼uX(δ) é transitiva. De fato, sejam x, x′, x′′ ∈ X tais que x ∼uX(δ) x ′ e x′ ∼uX(δ) x ′′. Então uX(x, x′), uX(x′, x′′) ≤ δ. Pela desigualdade triangular forte, temos que uX(x, x′′) ≤ max{uX(x, x′), uX(x′, x′′)} ≤ δ, de onde segue que x ∼uX(δ) x ′′. Agora, mostremos que DX = {DX(δ) | δ ≥ 0} é um dendrograma. (D1) Para δ = 0, temos que uX(x, x′) = 0 se, e somente se, x = x′, ∀x, x′ ∈ X. Logo, DX(0) = {{x}, x ∈ X}. Tomando δ0 = maxx 6=x′{uX(x, x′)}, temos que uX(x, x′) ≤ δ0, ∀x, x′ ∈ X. Logo, x ∼uX(δ0) x ′, ∀x, x′ ∈ X, isto é, DX(δ0) = {X}. (D2) Para δ1 ≤ δ2 e x, x′ ∈ X tais que x ∼uX(δ1) x ′, temos que uX(x, x′) ≤ δ1 ≤ δ2. Logo x ∼uX(δ2) x ′. (D3) Para cada δ ≥ 0 tal que DX(δ) 6= {X}, tomamos ε(δ) um número positivo qualquer tal que 0 < ε(δ) < m− δ, onde m = min x,x′∈X {uX(x, x′) | uX(x, x′) > δ}. Note que a �nitude de X garante a existência de ε(δ). Daí, para qualquer δ′ ∈ [δ, δ + ε(δ)], se x ∼uX(δ) x ′ então x ∼uX(δ′) x ′, ou seja, uX(x, x′) ≤ δ′ ≤ δ + ε(δ) < m. De fato, suponha, por absurdo, que uX(x, x′) > δ, então uX(x, x′) > m, o que é uma contradição. Para DX(δ) = {X}, a propriedade (D3) é satisfeita trivialmente, pois o dendro- grama se mantém inalterado para qualquer parâmetro su�cientemente grande. Em resumo, como ∼uX(δ) é relação de equivalência e DX é um dendrograma, segue que Υ está bem de�nida. Parte 2. Mostremos que Ψ ◦ Υ = id. Considere um conjunto ultramétrico (X, uX) e escreva δ0 = uX(x, x′), para x, x′ ∈ X dados. Sendo (X, u∗X) = Ψ ◦ Υ(X, uX), basta mostrar que uX(x, x′) = u∗X(x, x′). Para isso, se δ < δ0, x e x′ não se agrupam e para δ = δ0, se agrupam em um único cluster, pois uX(x, x′) = δ0. Assim, obtemos que δ0 = min{δ ≥ 0 | x ∼u∗X(δ) x ′}, ou seja, u∗X(x, x′) = δ0 = uX(x, x′). Finalmente, mostremos que Υ ◦Ψ = id. Para isso, considere DX um dendrograma com a relação de equivalência ∼DX(δ) e (X, uX) = Ψ(DX). Seja D∗X = Υ◦Ψ(DX), cuja relação de equivalência ∼D∗X proveniente da ultramétrica uX , para cada parâmetro δ ≥ 0, é dada por x ∼D∗X(δ) x ′ ⇔ uX(x, x′) ≤ δ. Assim, claramente as relações ∼DX(δ) e ∼D∗X(δ) são equivalentes, já que ambas são de�nidas a partir da mesma ultramétrica uX . Logo, os dendrogramas correspondentes são iguais, isto é, DX = D∗X . Dada a equivalência entre o conjunto dos dendrogramas e o das ultramétricas es- tabelecida no Teorema 2.3, agora podemos considerar os métodos de clustering hierár- quicos H com as ultramétricas induzidas no conjunto X baseado nas funções dissimi- laridades AX . Contudo, ultramétricas são casos particulares de funções dissimilaridades e assim podemos reinterpretar o método H como a função H : N → U 30 Axiomas do Valor e da Transformação onde u(x, x′) induzida por H é o parâmetro mínimo para o qual x e x′ serão agrupados por H, para x, x′ ∈ X. Dois métodos H1 e H2 são equivalentes se, e somente se, H1(N) = H2(N), para qualquer N ∈ N . (A1) Axioma do Valor. A imagem de H aplicada à network de dois nós satisfaz up,q(p, q) = max(α, β). (A2) Axioma da Transformação. Considere duas networks NX = (X,AX) e NY = (Y,AY ) e uma função redução de dissimilaridade φ : X → Y . Então para quais- quer x, x′ ∈ X, as imagens (X, uX) = H((X,AX)) e (Y, uY ) = H((Y,AY )) satis- fazem uX(x, x′) ≥ uY (φ(x), φ(x′)). Para o caso particular de networks simétricas, de�nimos o dendrograma single- linkage SLX através da relação de equivalência x ∼SLX(δ) x ′ ⇔ ũ∗X(x, x′) = ũ∗X(x′, x) ≤ δ. (2.4) Em outras palavras, para δ ≥ 0, o método single-linkage faz com que x e x′ façam parte do mesmo cluster se, e somente se, o custo mínimo de cadeia entre x e x′ não exceda δ. Pelo Teorema 2.3, existe uma equivalência entre o dendrograma proveniente da relação (2.4) e o conjunto ultramétrico, que denotamos por (X, uSLX ), em que uSLX (x, x′) = ũ∗X(x, x′) = ũ∗X(x′, x). (2.5) Para parâmetros 0 ≤ δ < mlc(X,AX) é impossível encontrar cadeias de in�uência mútua com custo máximo menor que δ entre quaisquer dois nós de X. De fato, su- ponhamos, por absurdo, que possamos ligar x a x′ em uma cadeia com custo máximo menor ou igual que δ. Da mesma forma, podemos ligar x′ a x. Assim, concatenando essas cadeias podemos formar um ciclo com custo menor ou igual a δ. Assim, temos a seguinte propriedade: (P1) Propriedade de In�uência. Para qualquer network NX = (X,AX), a imagem (X, uX) = H((X,AX)) correspondente à função do método de clustering hierár- quico é tal que a ultramétrica uX(x, x′) entre quaisquer dois nós x e x′ não pode ser menor que o custo mínimo de ciclo mlc(X,AX) da network, ou seja, uX(x, x′) ≥ mlc(X,AX), ∀x 6= x′. Uma segunda a�rmação intuitiva sobre in�uência em networks de tamanho arbi- trário é formalizada na forma do Axioma do Valor Estendido. Para isso, vamos introduzir o seguinte conceito. De�nição 2.4. De�na uma família canônica de networks assimétricas, representada na Figura 2.1, por −→ ∆n(α, β) := ({1, . . . , n}, An,α,β), (2.6) com n ∈ N e α, β > 0, e para i < j, o valor da dissimilaridade é An,α,β(i, j) = α e para i > j, a dissimilaridade é An,α,β(i, j) = β. Algumas equivalências 31 Na forma matricial, temos An,α,β =  0 α α α · · · α β 0 α α · · · α β β 0 α · · · α ... ... ... . . . ... ... β β β · · · 0 α β β β · · · β 0  1 2 . . . . . . n α α α β α α β β β Figura 2.1: Network canônica Agora, considere uma permutação Π = {π1, π2, . . . , πn} de {1, . . . , n} e a ação Π(A) de Π em uma matriz An×n, que de�nimos por (Π(A))i,j = Aπi,πj , ∀i, j. De�nindo a network −→ ∆n(α, β,Π) := ({1, . . . , n},Π(An,α,β)), podemos introduzir o conceito do Axioma do Valor Estendido. (A1′) Axioma do Valor Estendido. Seja −→ ∆n(α, β,Π) := ({1, . . . n},Π(An,α,β)). Então, para todos os índices n ∈ N, constantes α, β > 0, e permutações Π de {1, . . . , n}, a saída ({1, . . . , n}, un,α,β) = H( −→ ∆n(α, β,Π)) do método de clustering hierárquico H satisfaz un,α,β(i, j) = max(α, β), ∀i 6= j. Observe que o Axioma do Valor (A1) é um caso particular do Axioma do Valor Estendido (A1′) para n = 2. Além disso, note que mlc( −→ ∆n(α, β)) = max(α, β) pois, para formar um ciclo nessa network, é necessário atravessar pelo menos um link de custo α e pelo menos um link de custo β, como mostra a Figura 2.1. Como a permutação não altera o custo de um ciclo, então mlc( −→ ∆n(α, β,Π)) = mlc( −→ ∆n(α, β)) = max(α, β). (2.7) 2.1 Algumas equivalências Nos próximos teoremas, veremos que os Axiomas do Valor e do Valor Estendido são equivalentes quando o Axioma da Transformação for válido. Além disso, veremos que se o método H for admissível, então a Propriedade da In�uência se veri�ca. Para isso, começaremos demonstrando o seguinte lema. Lema 2.5. Seja N = (X,AX) uma network com n nós e δ > 0 uma constante. Suponha que x, x′ ∈ X são tais que o custo mínimo de cadeia associado satisfaz ũ∗X(x, x′) ≥ δ. (2.8) 32 Axiomas do Valor e da Transformação Então existe uma partição Pδ(x, x ′) = {Bδ(x), Bδ(x ′)} do conjunto X em clusters Bδ(x) 3 x e Bδ(x ′) 3 x′ tais que, para todos os pontos b ∈ Bδ(x) e b′ ∈ Bδ(x ′), AX(b, b′) ≥ δ. (2.9) Demonstração. Suponhamos, por absurdo, que não exista uma partição Pδ(x, x′) = {Bδ(x), Bδ(x ′)} com x ∈ Bδ(x) e x′ ∈ Bδ(x ′) satisfazendo (2.9), para quaisquer x, x′ ∈ X satisfazendo (2.8). Então existe pelo menos um par de nós x, x′ ∈ X satisfazendo (2.8) tal que para toda partição de X em dois clusters P = {B,B′} com x ∈ B e x′ ∈ B′, podemos encontrar pelo menos um par de elementos bP ∈ B e b′P ∈ B′ para os quais AX(bP , b ′ P ) < δ. (2.10) Consideremos, primeiramente, a partição P1 = {B1, B ′ 1}, onde B1 = {x} e B′1 = X − B1. Como (2.10) é verdadeiro para todas as partições e x é o único elemento de B1, então existe b′P1 ∈ B′1 tal que AX(x, b′P1 ) < δ. (2.11) Assim, a cadeia C(x, b′P1 ) = [x, b′P1 ] composta por esses dois nós tem custo menor que δ. Além disso, como ũ∗X(x, b′P1 ) representa o custo mínimo sobre todas as cadeias que ligam x a b′P1 , podemos concluir que ũ∗X(x, b′P1 ) ≤ AX(x, b′P1 ) < δ. A seguir, consideremos a partição P2 = {B2, B ′ 2}, onde B2 = {x, b′P1 } e B′2 = X−B2. Por (2.10), existe um nó b′P2 ∈ B′2 que satisfaz pelo menos uma das seguintes condições: AX(x, b′P2 ) < δ, (2.12) AX(b′P1 , b′P2 ) < δ. (2.13) Se (2.12) for verdadeira, a cadeia C(x, b′P2 ) = [x, b′P2 ] tem custo menor que δ. Se (2.13) for verdadeira, combinando com a a�rmação (2.11), temos que a cadeia C(x, b′P2 ) = [x, b′P1 , b′P2 ] tem custo menor que δ. Em ambos os casos, concluímos que existe uma cadeia C(x, b′P2 ) ligando x a b′P2 cujo custo é menor que δ. Consequente- mente, o custo mínimo de cadeia deve satisfazer ũ∗X(x, b′P2 ) < δ. Consideremos, agora, a partição P3 = {B3, B ′ 3} com B3 = {x, b′P1 , b′P2 } e B′3 = X − B3. Segue de (2.10) que existe um ponto b′P3 tal que pelo menos uma das dissi- milaridades AX(x, b′P3 ), AX(b′P1 , b′P3 ) ou AX(b′P2 , b′P3 ) é menor que δ, e assim como foi argumentado em (2.12) e (2.13), existe uma cadeia C(x, b′P3 ) ligando x a b′P3 cujo custo é menor que δ e daí segue que ũ∗X(x, b′P3 ) < δ. Essa construção pode ser repetida n−1 vezes, obtendo-se partições P1, P2, . . . , Pn−1 com os nós correspondentes b′P1 , b′P2 , . . . , b′Pn−1 tais que o custo mínimo de cadeia satisfaz ũ∗X(x, b′Pi) < δ, 1 ≤ i ≤ n− 1. (2.14) Observe que, por construção, os nós b′Pi são distintos entre si e também distintos de x. Como a network tem n nós, temos que x′ = b′Pk para algum k ∈ {1, . . . , n− 1}. Segue de (2.14) que ũ∗X(x, x′) < δ, o que é uma contradição, pois x e x′ satisfazem (2.8). Algumas equivalências 33 Teorema 2.6. Seja H um método de clustering hierárquico que satisfaz o Axioma da Transformação (A2). Então H satisfaz o Axioma do Valor (A1) se, e somente se, H satisfaz o Axioma do Valor Estendido (A1′). Demonstração. Suponhamos que H satisfaz o Axioma do Valor Estendido (A1′). Por- tanto H satisfaz o Axioma do Valor (A1), pois o Axioma (A1) é um caso particular do Axioma (A1′). Reciprocamente, considere H um método satisfazendo (A1) e (A2), −→ ∆n(α, β,Π) a network dada em (A1′) e ({1, . . . , n}, u) = H( −→ ∆n(α, β,Π)) sua imagem por H. Temos que mostrar que, para todo índice n ∈ N, constantes α, β > 0, permutações Π de {1, . . . , n} e pontos i 6= j, un,α,β(i, j) = max(α, β). Assim, mostraremos un,α,β(i, j) ≤ max(α, β), (2.15) un,α,β(i, j) ≥ max(α, β). (2.16) Para (2.15), de�na uma network de dois nós simétrica −→ ∆2(max(α, β),max(α, β)) = ({p, q}, Ap,q), onde Ap,q(p, q) = Ap,q(q, p) = max(α, β) e denote por ({p, q}, up,q) = H( −→ ∆2(max(α, β),max(α, β))). Como H satisfaz (A1), up,q(p, q) = max(max(α, β),max(α, β)) = max(α, β). (2.17) Considere, agora, a função φi,j : {p, q} → {1, . . . , n}, onde φi,j(p) = i e φi,j(q) = j. Como as dissimilaridades em −→ ∆n(α, β,Π) são α ou β e as dissimilaridades em −→ ∆2(max(α, β),max(α, β)) são max(α, β), segue que a função φi,j é uma redução de dissimilaridade, independentemente dos valores i, j. Como o método H satisfaz (A2), temos que up,q(p, q) ≥ un,α,β(φi,j(p), φi,j(q)) = un,α,β(i, j). (2.18) Substituindo (2.17) em (2.18), segue a desigualdade (2.15). Para mostrar a desigualdade (2.16), considere dois nós arbitrários distintos i, j ∈ {1, . . . , n} da network −→ ∆n(α, β,Π). Denote por C(i, j) e C(j, i) duas cadeias cujo custo é mínimo, ou seja, são as cadeias que realizam o custo mínimo de cadeias ũ∗n,α,β(i, j) e ũ∗n,α,β(j, i), respectivamente. Observe que deve ser verdadeira pelo menos uma das desigualdades: ũ∗n,α,β(i, j) ≥ max(α, β), (2.19) ũ∗n,α,β(j, i) ≥ max(α, β). (2.20) De fato, se ambas forem falsas, a concatenação de C(i, j) e C(j, i) formaria um ciclo C(i, i) = C(i, j) ] C(j, i) de custo estritamente menor que max(α, β), o que não pode ocorrer, pois mlc( −→ ∆n(α, β,Π)) = max(α, β) por (2.7). Sem perda de generalidade, assuma que vale (2.19) e considere δ = max(α, β). Pelo Lema 2.5, podemos garantir que existe uma partição do conjunto de nós {1, . . . , n} em dois blocos Bδ(i) e Bδ(j), com i ∈ Bδ(i) e j ∈ Bδ(j) tais que para todo b ∈ Bδ(i) e b′ ∈ Bδ(j), temos Π(An,α,β)(b, b′) ≥ δ = max(α, β). (2.21) De�na uma network de dois nós −→ ∆2(max(α, β),min(α, β)) = ({r, s}, Ar,s), em que Ar,s(r, s) = max(α, β) e Ar,s(s, r) = min(α, β) e escreva ({r, s}, ur,s) = H( −→ ∆2(max(α, β),min(α, β))). 34 Axiomas do Valor e da Transformação Como H satisfaz (A1), temos ur,s(r, s) = max(max(α, β),min(α, β)) = max(α, β). (2.22) Considere a função φ′i,j : {1, . . . , n} → {r, s} tal que φ′i,j(b) = r e φ′i,j(b ′) = s, para todo b ∈ Bδ(i) e b′ ∈ Bδ(j). A função φ′i,j é redução de dissimilaridade, pois Π(An,α,β)(k, l) ≥ Ar,s(φ ′ i,j(k), φ′i,j(l)), ∀k, l ∈ {1, . . . , n}. (2.23) Com efeito, considere três possíveis casos. � Se k e l pertencem a um mesmo bloco, isto é, k, l ∈ Bδ(i) ou k, l ∈ Bδ(j), então φ′i,j(k) = φ′i,j(l) e Ar,s(φ′i,j(k), φ′i,j(l)) = 0, e portanto não pode ser maior que Π(An,α,β)(k, l). � Se k ∈ Bδ(j) e l ∈ Bδ(i), segue que Ar,s(φ′i,j(k), φ′i,j(l)) = Ar,s(s, r) = min(α, β). Logo, não pode exceder Π(An,α,β)(k, l), que é igual a α ou β. � Se k ∈ Bδ(i) e l ∈ Bδ(j), temos Ar,s(φ′i,j(k), φ′i,j(l)) = Ar,s(r, s) = max(α, β), mas por (2.21), temos que Π(An,α,β)(k, l) ≥ δ = max(α, β). Como H satisfaz (A2), segue que un,α,β(i, j) ≥ ur,s(φ ′ i,j(i), φ ′ i,j(j)) = ur,s(r, s). (2.24) Substituindo (2.22) em (2.24), obtemos a desigualdade (2.16), o que completa a prova. Lema 2.7. Sejam N = (X,AX) uma network arbitrária com n nós e −→ ∆n(α, β) = ({1, . . . , n}, An,α,β) a network canônica com 0 < α ≤ sep(X,AX) e β = mlc(X,AX). Então existe uma bijeção φ : X → {1, . . . , n} tal que AX(x, x′) ≥ An,α,β(φ(x), φ(x′)), ∀x, x′ ∈ X. (2.25) Demonstração. Para construir a bijeção φ, consideremos a função P : X → P(X) tal que P (x) = {x′ ∈ X | x′ 6= x,AX(x′, x) < β}. É importante observar que existe x ∈ X tal que P (x) = ∅. De fato, suponhamos, por absurdo, que tal fato não ocorra. Daí, considere xn ∈ X um nó qualquer e construa a cadeia [x0, x1, . . . , xn] onde o i-ésimo elemento da cadeia, xi−1, está na P -imagem de xi. Da de�nição da função P , segue que todas as dissimilaridades dessa cadeia satisfazem AX(xi−1, xi) < β = mlc(X,AX). Mas como a cadeia [x0, x1, . . . , xn] contém n + 1 elementos, pelo menos um nó deve ser repetido, e assim encontramos um ciclo para o qual todas as dissimilaridades estão limitadas por β = mlc(X,AX), o que é impossível, pois contradiz a de�nição de mlc. Portanto, existe um nó xi1 para o qual P (xi1) = ∅. De�nimos φ(xi1) = 1. Tomemos agora xi2 6= xi1 cuja P -imagem é {xi1} ou ∅, ou seja, P (xi2) ⊆ {xi1}. Tal nó existe, pois caso contrário podemos tomar xn−1 ∈ X − {xi1} e construir a cadeia [x0, x1, . . . , xn−1], em que xi−1 ∈ P (xi) − {xi1}, isto é, xi−1 está na P -imagem de xi e xi−1 6= xi1 . Como a cadeia [x0, x1, . . . , xn−1] contém n elementos do conjunto X−{xi1} de cardinalidade n−1, pelo menos um dos nós se repete. Assim, encontramos um ciclo Algumas equivalências 35 onde todas as dissimilaridades entre nós consecutivos satisfazem AX(xi−1, xi) < β = mlc(X,AX), contradizendo a de�nição de mlc. De�nimos φ(xi2) = 2. Repetindo esse processo, no k-ésimo passo temos φ(xik) = k para algum nó xik /∈ {xi1 , xi2 , . . . , xik−1 }, cuja P -imagem é um subconjunto dos nós já escolhidos, ou seja, P (xik) ⊆ {xi1 , xi2 , . . . , xik−1 }. Como todos os nós xik são distintos, para todo k, a função φ dada por φ(xik) = k é bijetiva. Por construção, para todo l > k, tem-se xil /∈ P (xik), e da de�nição de P , vale AX(xil , xik) ≥ β, ∀l > k. Além disso, da de�nição da dissimilaridade An,α,β, segue que An,α,β(φ(xil), φ(xik)) = An,α,β(l, k) = β, ∀l > k. Assim, concluímos que (2.25) é verdadeiro para todos os nós em que φ(x) > φ(x′). Quando φ(x) < φ(x′), temos que An,α,β(φ(x), φ(x′)) = α ≤ sep(X,AX) ≤ mlc(X,AX) = β. Portanto AX(x, x′) ≥ An,α,β(φ(x), φ(x′)), ∀x, x′ ∈ X. Teorema 2.8. Considere H um método de clustering hierárquico que satisfaz (A1′) e (A2). Então H satisfaz (P1). Demonstração. Seja N = (X,AX) uma network arbitrária tal que X = {x1, . . . , xn} e denote por (X, uX) = H(X,AX). Assim, basta provar que uX(x, x′) ≥ mlc(X,AX), para todo x 6= x′. Considere a network canônica −→ ∆n(α, β) = ({1, . . . , n}, An,α,β) com β = mlc(X,AX) e 0 < α ≤ sep(X,AX). Assim, temos que 0 < α ≤ sep(X,AX) ≤ mlc(X,AX) = β. Denote por ({1, . . . , n}, uα,β) = H( −→ ∆n(α, β)). Como H satisfaz (A1′), então para todos os índices i, j ∈ {1, . . . , n} tal que i 6= j, temos uα,β(i, j) = max(α, β) = β = mlc(X,AX). (2.26) Além disso, considere a função redução de dissimilaridade bijetiva do Lema 2.7 e note que, como o método H satisfaz (A2), para todo x, x′ ∈ X vale uX(x, x′) ≥ uα,β(φ(x), φ(x′)). Como a igualdade em (2.26) é verdadeira para todo i 6= j e como φ é bijetiva, temos que uX(x, x′) ≥ β = mlc(X,AX), para x 6= x′ em X, e portanto H satisfaz (P1). Corolário 2.9. Considere H um método de clustering hierárquico que satisfaz (A1) e (A2), então H satisfaz (P1). Demonstração. Pelo Teorema 2.6, temos que se H satisfaz (A2), então (A1) e (A1′) são equivalentes, e portanto, pelo Teorema 2.8, H satisfaz (P1). 3 Clustering Recíproco e Não-Recíproco Neste capítulo, iremos construir dois métodos de clustering hierárquico que satis- fazem (A1) e (A2), chamados de clustering recíproco e não-recíproco. Seja NX = (X,AX) ∈ N uma network qualquer. Considere a dissimilaridade simétrica ĀX(x, x′) := max{AX(x, x′), AX(x′, x)}, ∀x, x′ ∈ X. (3.1) De�nição 3.1. De�nimos o método de clustering recíproco HR por (X, uRX) = HR(X,AX), onde a ultramétrica uRX(x, x′) entre os nós x e x′ é dada por uRX(x, x′) := min C(x,x′) max i|xi∈C(x,x′) ĀX(xi, xi+1). (3.2) Lembremos que, pelo Teorema 2.3, existe uma equivalência entre os dendrogramas e as ultramétricas. Sendo RX o dendrograma produzido pelo clustering recíproco, podemos escrever as classes de equivalência do clustering recíproco como x ∼RX(δ) x ′ ⇔ min C(x,x′) max i|xi∈C(x,x′) ĀX(xi, xi+1) ≤ δ. Proposição 3.2. O método de clustering recíproco HR é válido e admissível. Demonstração. Para provar que HR é válido e admissível, temos que mostrar que uRX é de fato uma ultramétrica e que valem os axiomas (A1) e (A2). Primeiro, mostremos que uRX é uma ultramétrica. As propriedades (i) e (ii) da De- �nição 2.2 seguem direto da de�nição de uRX . Para veri�car a desigualdade triangular forte, considere C∗(x, x′) e C∗(x′, x′′) duas cadeias que atingem o mínimo em (3.2) para uRX(x, x′) e uRX(x′, x′′), respectivamente. Assim, o custo máximo na cadeia concatenada C(x, x′′) = C∗(x, x′) ]C∗(x′, x′′) não excede o custo máximo em cada cadeia individu- almente. Portanto, embora o custo máximo possa ser menor em cadeias diferentes, a cadeia C(x, x′′) é su�ciente para que uRX(x, x′′) ≤ max(uRX(x, x′), uRX(x′, x′′)). Agora, provemos que HR satisfaz os Axiomas (A1) e (A2). Dada uma network de dois nós −→ ∆2(α, β), denote ({p, q}, uRp,q) = HR( −→ ∆2(α, β)). Como toda cadeia possível de p para q deve conter p e q como nós consecutivos, temos que uRp,q(p, q) = max(Ap,q(p, q), Ap,q(q, p)) = max(α, β), 37 38 Clustering Recíproco e Não-Recíproco provando, assim, o Axioma (A1). A seguir, considere duas networks (X,AX) e (Y,AY ) quaisquer e uma função re- dução de dissimilaridade φ : X → Y . Sejam (X, uRX) = HR(X,AX) e (Y, uRY ) = HR(Y,AY ). Denote por C∗X(x, x′) = [x = x0, x1, . . . , xl = x′] a cadeia que atinge o mínimo em (3.2), com x, x′ ∈ X arbitrários e considere CY (φ(x), φ(x′)) = [φ(x) = φ(x0), φ(x1), . . . , φ(xl) = φ(x′)], uma cadeia em Y . Como a função φ é uma redução de dissimilaridade, temos que AY (φ(xi), φ(xi+1)) ≤ AX(xi, xi+1) e AY (φ(xi+1), φ(xi)) ≤ AX(xi+1, xi), de onde obtemos max φ(xi)∈CY (φ(x),φ(x′)) ĀY (φ(xi), φ(xi+1)) ≤ uRX(x, x′). Além disso, note que CY (φ(x), φ(x′)) é uma cadeia particular que liga φ(x) a φ(x′) e que a ultramétrica recíproca faz uso da cadeia de custo mínimo entre todas as cadeias que ligam φ(x) a φ(x′). Portanto, uRY (φ(x), φ(x′)) ≤ max φ(xi)∈CY (φ(x),φ(x′)) ĀY (φ(xi), φ(xi+1)) ≤ uRX(x, x′), o que conclui a prova. De�nição 3.3. Seja NX = (X,AX) ∈ N uma network qualquer. De�nimos o método de clustering não-recíproco HNR por (X, uNRX ) = HNR(X,AX), onde a ultramétrica uNRX (x, x′) entre os nós x e x′ é dada por uNRX (x, x′) := max(ũ∗X(x, x′), ũ∗X(x′, x)), (3.3) lembrando que ũ∗X(x, x′) = minC(x,x′) maxi|xi∈C(x,x′)AX(xi, xi+1). Agora, provaremos uma proposição análoga à Proposição 3.2 para clustering não- recíproco. Proposição 3.4. O método de clustering não-recíproco HNR é válido e admissível. Demonstração. Como na proposição anterior, iniciamos mostrando a desigualdade tri- angular forte para uNRX . Para isso, considere as cadeias C∗(x, x′) e C∗(x′, x′′) que atin- gem o mínimo em ũ∗X(x, x′) e ũ∗X(x′, x′′), bem como as cadeias C∗(x′′, x′) e C∗(x′, x) que atingem o mínimo em ũ∗X(x′′, x′) e ũ∗X(x′, x), respectivamente. As cadeias concatenadas C(x, x′′) = C∗(x, x′)]C∗(x′, x′′) e C∗(x′′, x) = C∗(x′′, x′)] C∗(x′, x) não podem exceder o custo máximo em cada cadeia individualmente, e por- tanto uNRX (x, x′′) = max(ũ∗X(x, x′′), ũ∗X(x′′, x)) ≤ max(ũ∗X(x, x′), ũ∗X(x′, x′′), ũ∗X(x′′, x′), ũ∗X(x′, x)) = max(uNRX (x, x′), uNRX (x′, x′′)). Agora mostremos que valem os Axiomas (A1) e (A2). Considere a network de dois nós −→ ∆2(α, β) e denote ({p, q}, uNRp,q ) = HNR( −→ ∆2(α, β)). Note que ũ∗p,q(p, q) = α e ũ∗p,q(q, p) = β. Assim, temos uNRp,q (x, x′) = max(ũ∗p,q(p, q), ũ ∗ p,q(q, p)) = max(α, β). Clustering Recíproco e Não-Recíproco 39 Finalmente, considere duas networks (X,AX) e (Y,AY ) quaisquer e uma função redução de dissimilaridade φ : X → Y . Sejam (X, uNRX ) = HNR(X,AX) e (Y, uNRY ) = HNR(Y,AY ). Denote por C∗X(x, x′) = [x = x0, x1, . . . , xl = x′] a cadeia que atinge o mínimo em ũ∗X(x, x′), com x, x′ ∈ X arbitrários e seja CY (φ(x), φ(x′)) = [φ(x) = φ(x0), φ(x1), . . . , φ(xl) = φ(x′)] uma cadeia de Y que liga φ(x) a φ(x′). Como a função φ é redução de dissimilaridade, temos que AY (φ(xi), φ(xi+1)) ≤ AX(xi, xi+1). Consequentemente, max φ(xi)∈CY (φ(x),φ(x′)) AY (φ(xi), φ(xi+1)) ≤ max i|xi∈C∗(x,x′) AX(xi, xi+1) = ũ∗X(x, x′). Como ũ∗Y (φ(x), φ(x′)) é o custo mínimo de cadeia, segue que ũ∗Y (φ(x), φ(x′)) ≤ max φ(xi)∈CY (φ(x),φ(x′)) AY (φ(xi), φ(xi+1)) ≤ max i|xi∈C∗(x,x′) AX(xi, xi+1) = ũ∗X(x, x′). (3.4) Como x e x′ são arbitrários, então (3.4) é válido para o par (x′, x), ou seja, ũ∗Y (φ(x′), φ(x)) ≤ ũ∗X(x′, x). (3.5) Por (3.2) e (3.4), segue que uNRY (φ(x), φ(x′)) = max(ũ∗Y (φ(x), φ(x′)), ũ∗Y (φ(x′), φ(x))) ≤ max(ũ∗X(x, x′), ũ∗X(x′, x)) = uNRX (x, x′), o que conclui a demonstração. Observação 3.5. Podemos observar nas Figuras 3.1 e 3.2 os métodos de Clustering recíproco e não-recíproco, respectivamente. No método recíproco, uRX(x, x′) ≤ δ se existirem cadeias (recíprocas) cujas dissimilaridades máximas são menores ou iguais a δ em ambas as direções. No método não-recíproco, uNRX (x, x′) ≤ δ se x e x′ se juntam em ambas as direções com cadeias (possivelmente diferentes) cujas dissimilaridades máximas são menores ou iguais a δ. Portanto, clustering recíproco é um caso particular de clustering não-recíproco, e assim, para quaisquer x, x′ ∈ X, uNRX (x, x′) ≤ uRX(x, x′). (3.6) x x1 . . . . . . xl−1 x′ AX(x, x1) AX(x1, x2) AX(x1, x) AX(x2, x1) AX(xl−2, xl−1) AX(xl−1, xl−2) AX(xl−1, x′) AX(x′, xl−1) Figura 3.1: Representação de um clustering recíproco O próximo teorema nos diz que, dado um método de clustering hierárquico admis- sível qualquer, tem-se que a ultramétrica proveniente desse método sempre permanece entre as ultramétricas uNRX e uRX . 40 Clustering Recíproco e Não-Recíproco x x1 x′l′−1 . . . . . . . . . . . . xl−1 x′1 x′ AX(x, x1) AX(x1, x2) AX(x′l′−1, x) AX(xl−2, xl−1) AX(x′l′−2, x ′ l′−1) AX(xl−1, x′) AX(x′1, x ′ 2) AX(x′, x′1) Figura 3.2: Representação de um clustering não-recíproco Teorema 3.6. Seja H um método de clustering hierárquico que satisfaz os axiomas (A1) e (A2). Dada uma network arbitrária N = (X,AX), se (X, uX) = H(N) então uNRX (x, x′) ≤ uX(x, x′) ≤ uRX(x, x′), para quaisquer nós x, x′ ∈ X. Demonstração. Faremos a demonstração em duas partes, uma para cada desigualdade. Parte 1. uNRX (x, x′) ≤ uX(x, x′), ∀x, x′ ∈ X. Seja ∼NRX(δ) a relação de equivalência do clustering não-recíproco, dada por x ∼NRX(δ) x ′ ⇔ uNRX (x, x′) ≤ δ. (3.7) Agora, considere o conjunto Z = X mod ∼NRX(δ) das classes de equivalência corres- pondentes à relação ∼NRX(δ) e a função φδ : X → Z dada por φδ(x) = z, em que z é a classe de equivalência correspondente a x. Logo, se φδ(x) = φδ(x ′), então x ∼NRX(δ) x ′ e, por (3.7), uNRX (x, x′) ≤ δ. Assim, podemos escrever φδ(x) = φδ(x ′)⇔ uNRX (x, x′) ≤ δ. (3.8) De�nimos a network NZ = (Z,AZ), em que AZ(z, z′) := min x∈φ−1 δ (z) x′∈φ−1 δ (z′) AX(x, x′). (3.9) Observe que, por construção, φδ é redução de dissimilaridade, isto é, AX(x, x′) ≥ AZ(φδ(x), φδ(x ′)). (3.10) De fato, se φδ(x) = φδ(x ′), temos que AZ(φδ(x), φδ(x ′)) = 0 e se φδ(x) 6= φδ(x ′), temos que minx∈φ−1 δ (z),x′∈φ−1 δ (z′)AX(x, x′) ≤ AX(x, x′). Logo, em ambos os casos, vale (3.10). Seja H um método de clustering hierárquico que satisfaz os axiomas (A1) e (A2) e denote por (Z, uZ) = H(NZ). Notemos que de acordo com a de�nição (3.9), se z 6= z′ então AZ(z, z′) > δ ou AZ(z′, z) > δ. De fato, se ambos AZ(z, z′) ≤ δ e AZ(z′, z) ≤ δ, poderemos construir cadeias C(x, x′) e C(x′, x) com custo máximo menor que δ para quaisquer x ∈ φ−1δ (z) e x′ ∈ φ−1δ (z′). Clustering Recíproco e Não-Recíproco 41 Para a cadeia C(x, x′), denote por xo ∈ φ−1δ (z) e x′i ∈ φ−1δ (z′) os nós que atingem o mínimo em (3.9), ou seja, AZ(z, z′) = AX(xo, x ′ i). Como x e xo estão na mesma classe de equivalência, existe uma cadeia C(x, xo) de custo máximo menor ou igual a δ. Da mesma forma, como x′ e x′i estão na mesma classe de equivalência, existe uma cadeia C(x′i, x ′) de custo máximo menor ou igual a δ. Consequentemente, a cadeia concatenada C(x, x′) = C(x, xo) ] [xo, x ′ i] ] C(x′i, x ′) (3.11) tem custo máximo menor ou igual a δ. Analogamente, podemos construir uma cadeia C(x′, x) cujo custo máximo é menor ou igual a δ. Contudo, a existência dessas duas cadeias implicaria que x e x′ estão no mesmo cluster no parâmetro δ, contrariando o fato de z ser diferente de z′. A�rmação: mlc(NZ) > δ. Suponhamos que mlc(NZ) ≤ δ. Denote por [z, z′, . . . , z(l), z] um ciclo cujo custo seja menor ou igual a δ. Para quaisquer x ∈ φ−1δ (z) e x′ ∈ φ−1δ (z′), podemos ligar x a x′ usando a cadeia concatenada (3.11). Para ligar x′ a x, denote por x(k)o e x(k+1) i os pontos para os quais AZ(z(k), z(k+1)) = AX(x (k) o , x (k+1) i ). Assim, podemos ligar x′o e x (l) o com a seguinte cadeia concatenada, C(x′o, x (l) o ) = l−i⊎ k=1 [x(k)o , x (k+1) i ] ] C(x (k+1) i , x(k+1) o ). O custo máximo dessa cadeia é menor ou igual a δ, pois o custo máximo na cadeia C(x (k+1) i , x (k+1) o ) é menor ou igual a δ, já que ambos os nós pertencem à mesma classe z(k+1), e AX(x (k) o , x (k+1) i ) ≤ δ, para todo k. Assim, podemos ligar x′ a x com a cadeia concatenada C(x′, x) = C(x′, x′o) ] C(x′o, x (l) o ) ] [x(l)o , xi] ] C(xi, x), cujo custo máximo é menor ou igual a δ. Usando as cadeias C(x, x′) e C(x′, x) construídas acima, podemos concluir que uNRX (x, x′) ≤ δ, o que contradiz o fato de que x e x′ pertencem a diferentes classes de equivalência. Consequentemente, mlc(NZ) > δ, como a�rmado. Segue da Propriedade da In�uência (P1) (e da a�rmação acima) que uZ(z, z′) > δ, ∀z, z′ ∈ Z, z 6= z′. Além disso, por (3.8), (3.10) e (A2), tem-se uX(x, x′) ≥ uZ(z, z′) > δ, ou seja, quando x e x′ são aplicadas em classes de equivalência distintas, temos que uX(x, x′) > δ, ou ainda, uNRX (x, x′) > δ. Consequentemente, podemos a�rmar que uNRX (x, x′) > δ implica uX(x, x′) > δ, ou seja, {(x, x′) | uNRX (x, x′) > δ} ⊆ {(x, x′) | uX(x, x′) > δ}. (3.12) Como (3.12) vale para δ > 0 arbitrário, segue que uNRX (x, x′) ≤ uX(x, x′), ∀x, x′ ∈ X. Parte 2. uX(x, x′) ≤ uRX(x, x′), ∀x, x′ ∈ X. 42 Clustering Recíproco e Não-Recíproco Sejam x e x′ arbitrários com uRX(x, x′) = δ. Seja C(x, x′) = [x = x0, . . . , xl = x′] uma cadeia que atinge o mínimo em (3.2). Assim, podemos escrever δ = uRX(x, x′) = max i {max(AX(xi, xi+1), AX(xi+1, xi))}. (3.13) Considere a network simétrica −→ ∆2(δ, δ) = ({p, q}, Ap,q) e sua imagemH( −→ ∆2(δ, δ)) = ({p, q}, up,q). Por (A1), temos que up,q(p, q) = max(δ, δ) = δ. Sejam φi : {p, q} → X as transformações dadas por φi(p) = xi, φi(q) = xi+1, para xi e xi+1 nós consecutivos da cadeia C(x, x′). Como AX(xi, xi+1) ≤ δ e AX(xi+1, xi) ≤ δ por (3.13), obtemos AX(φi(p), φi(q)) ≤ Ap,q(p, q) = δ, AX(φi(q), φi(p)) ≤ Ap,q(q, p) = δ, o que implica que φi é redução de dissimilaridade e daí segue de (A2) que uX(φi(p), φi(q)) ≤ up,q(p, q) = δ. (3.14) Como (3.14) é válida para todo i e φi(p) = xi, φi(q) = xi+1, concluímos que uX(xi, xi+1) ≤ δ, ∀i. Pela desigualdade triangular forte, temos que uX(x, x′) ≤ max i {uX(xi, xi+1)} ≤ δ = uRX(x, x′), ∀x, x′ ∈ X. Networks simétricas Finalizamos este capítulo trabalhando apenas com networks simétricas, ou seja, vamos considerar N = (X,AX) tal que AX(x, x′) = AX(x′, x), para todo x, x′ ∈ X. Notemos que nesse caso os clusterings recíproco e não-recíproco são equivalentes. De fato, dados x, x′ ∈ X arbitrários, ĀX(xi, xi+1) = AX(xi, xi+1) = AX(xi+1, xi), assim temos que a de�nição da ultramétrica recíproca se reduz a uRX(x, x′) = min C(x,x′) max i|xi∈C(x,x′) ĀX(xi, xi+1) = min C(x,x′) max i|xi∈C(x,x′) AX(xi, xi+1) = ũ∗X(x, x′). Além disso, note que os custos de uma cadeia C(x, x′) = [x = x0, . . . , xl = x′] e de sua recíproca C(x′, x) = [x′ = xl, . . . , x0 = x] são os mesmos, o que implica que ũ∗X(x, x′) = ũ∗X(x′, x) e, portanto, pela de�nição de ultramétrica não-recíproca, temos uNRX (x, x′) = ũ∗X(x, x′) = uRX(x, x′). (3.15) Lembremos de (2.5) que a ultramétrica single-linkage para networks simétricas satisfaz uSLX (x, x′) = ũ∗X(x, x′), e portanto, por (3.15), temos que uSLX (x, x′) = uRX(x, x′) = uNRX (x, x′). (3.16) Assim, pelo Teorema 3.6, concluímos que uSLX é a única ultramétrica possível para networks simétricas. Mais precisamente: Corolário 3.7. Sejam H : N → U um método de clustering hierárquico qualquer, N = (X,AX) uma network simétrica e HSL o método single-linkage. Se H satisfaz os axiomas (A1) e (A2), então H = HSL. Demonstração. Como H satisfaz as hipóteses do Teorema 3.6, temos que uNRX (x, x′) ≤ uX(x, x′) ≤ uRX(x, x′), para todo x, x′ ∈ X. Mas, por (3.16), uSLX (x, x′) = uRX(x, x′) = uNRX (x, x′), o que implica que uSLX (x, x′) = uX(x, x′), para todo x, x′ ∈ X, ou seja, H ≡ HSL. 4 Algoritmos Neste capítulo, interpretaremos as de�nições de ultramétrica recíproca e não-recí- proca por meio de matrizes de modo a resolvermos problemas computacionalmente. Para isso, introduziremos primeiramente o conceito de dioide, necessário para tal in- terpretação. Mais detalhes sobre Álgebra Dioide podem ser encontrados em [5]. 4.1 Álgebra Dioide Seja E um conjunto não-vazio qualquer e ⊕ uma operação binária em E. Dizemos que (E,⊕) é um monoide se ⊕ for associativa e tiver um elemento neutro. Se ⊕ for comutativa, dizemos que o monoide (E,⊕) é comutativo. Um elemento a ∈ E é idempotente se a ⊕ a = a. Além disso, ⊕ é idempotente se a⊕ a = a,∀a ∈ E. De�nição 4.1. Dizemos que (E,⊕,⊗) é um dioide se: (i) (E,⊕) é um monoide comutativo com elemento neutro ε; (ii) (E,⊗) é um monoide com elemento neutro e; (iii) A relação a ≤ b ⇔ ∃c : b = a ⊕ c é uma relação de ordem, isto é, se a ≤ b e b ≤ a então a = b; (iv) ε é absorvente para ⊗, ou seja, para qualquer a ∈ E, a⊗ ε = ε⊗ a = ε; (v) ⊗ é distributiva com respeito a ⊕. Considere (E,⊕,⊗) um dioide qualquer. Como E é um conjunto ordenado dado na De�nição 4.1, item (iii), podemos estabelecer uma topologia para E induzida por essa relação de ordem. Para mais detalhes, veja em [5, Capítulo 3, Seção 2]. Isto posto, podemos falar sobre limite e continuidade. De�nição 4.2. Sejam (E,⊕,⊗) um dioide e a ∈ E qualquer. O quasi-inverso de a, denotado por a∗, é o limite, quando existir, da sequência a∗ := lim k→∞ e⊕ a⊕ a(2) ⊕ · · · ⊕ a(k), onde e é o elemento neutro de ⊗ e a(k) = a⊗ · · · ⊗ a︸ ︷︷ ︸ k vezes . 43 44 Algoritmos Um elemento a ∈ E é p-estável se e⊕ a⊕ a(2) ⊕ · · · ⊕ a(p) = e⊕ a⊕ a(2) ⊕ · · · ⊕ a(p) ⊕ a(p+1) = · · · onde e é o elemento neutro de ⊗. A De�nição 4.2 pode ser interpretada utilizando-se matrizes, como a seguir. A quasi-inversa de uma matriz A, denotada por A∗, é o limite, quando existir, da sequência de matrizes A∗ := lim k→∞ I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(k), (4.1) onde I é a matriz identidade na álgebra dioide, ou seja, I tem o elemento neutro e na diagonal e o elemento neutro ε fora da diagonal. Antes de prosseguirmos, faremos brevemente menção a alguns conceitos e resulta- dos, obtidos com mais detalhes em [5, Capítulo 4, Seção 3.2]. Proposição 4.3. Se ⊕ é idempotente, então (I ⊕ A)(k) = I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(k). Demonstração. Temos que (I ⊕ A)(k) = (I ⊕ A)⊗ · · · ⊗ (I ⊕ A)︸ ︷︷ ︸ k vezes = I ⊕ k⊕ r=1 Cr kA (r), (4.2) onde Cr k = k! r!(k − r)! . Como ⊕ é idempotente e Cr k ∈ Z, temos que Cr kA (r) = A(r). Substituindo em (4.2), concluímos a prova. Sejam (E,⊕,⊗) um dioide, A uma matriz n× n com entradas em E e G(A) o seu grafo correspondente, ou seja, tendo A como matriz de adjacência. Um circuito de G(A) é um circuito baseado quando um de seus vértices é tratado como origem. Neste caso, o peso do circuito baseado γ = {i1i2, i2i3, . . . , iki1} de origem i1 é w(γ) = ai1i2 ⊗ ai2i3 ⊗ · · · ⊗ aiki1 . Dizemos que um grafo G(A) não tem um circuito p-absorvente se o peso de cada circuito baseado é um elemento p-estável de E. Com isso, podemos apresentar (ver [5, Capítulo 3, Seção 3.3, Teorema 1]): Teorema 4.4. Se G(A) não tem um circuito 0-absorvente, então A∗ = lim k→+∞ I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(k) = I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(n−1) = I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(n) = · · · Além disso, A∗ satisfaz a equação A∗ = I ⊕ A⊗ A∗ = I ⊕ A∗ ⊗ A. 4.2 Algoritmo para ultramétricas Dada uma network N = (X,AX) tal que |X| = n, interpretaremos AX como uma matriz de dissimilaridades n× n. Também, podemos denotar a ultramétrica uX como uma matriz simétrica n× n e reinterpretar (3.1) como ĀX := max(AX , A T X), Algoritmo para ultramétricas 45 em que ATX é a matriz transposta de AX e max(−,−) é tomado ponto-a-ponto, ou seja, resulta em uma matriz de máximos. Assim, podemos interpretar as de�nições de ultramétricas recíproca e não-recíproca como matrizes e então resolver problemas computacionalmente. Para isso, usaremos os conceitos vistos na Seção 4.1 para a álgebra dioide U := (R+ ∪ {+∞},min,max), conforme em [1]. Nesse caso, a soma ⊕ é a operação mínimo e o produto ⊗ é a operação máximo. Mais precisamente, dados a, b ∈ R+ ∪ {+∞}, então a⊕ b = min(a, b) e a⊗ b = max(a, b). Denotaremos o conjunto {1, 2, . . . , n} por [1, n]. Na álgebra dioide, o produto A ⊗ B de duas matrizes A e B é dado pela matriz com entradas [A⊗B]ij := n⊕ k=1 (Aik ⊗Bkj) = min k∈[1,n] max(Aik, Bkj). (4.3) Para inteiros k ≥ 2, a k-ésima potência de uma matriz de dissimilaridades AX relacionada à matriz ultramétrica uX é dada por A (k) X := AX ⊗ A(k−1) X , com A (1) X := AX . Note que os elementos do produto dioide u(2)X são dados por [u (2) X ]ij = min k∈[1,n] max([uX ]ik, [uX ]kj). Como uX satisfaz a desigualdade triangular forte, então [uX ]ij ≤ max([uX ]ik, [uX ]kj), para todo k ∈ [1, n]. Em particular, para k = j, temos max([uX ]ij, [uX ]jj) = max([uX ]ij, 0) = [uX ]ij. Consequentemente, [u (2) X ]ij = [uX ]ij, para todo i, j ∈ [1, n], e portanto u(2)X = uX , ou seja, uX é idempotente. Na de�nição de quasi-inversa de uma matriz A, dada em (4.1), A∗ := lim k→∞ I ⊕ A⊕ A(2) ⊕ · · · ⊕ A(k). Porém, para o dioide U := (R+ ∪ {+∞},min,max), a matriz identidade I tem 0 na diagonal e +∞ fora da diagonal. A utilidade de uma quasi-inversa para o nosso estudo reside no fato de que, dada uma matriz dissimilaridade AX , tem-se [A∗X ]ij = min C(xi,xj) max k|xk∈C(xi,xj) AX(xk, xk+1), (4.4) onde A∗X é a quasi-inversa de AX . Em outras palavras, os elementos da quasi-inversa A∗X correspondem aos custos mínimos de cadeia da network associada (X,AX). Para mais detalhes, veja [5, Capítulo 6, Seção 6.1]. Assim, para calcular computacionalmente as ultramétricas recíproca e não-recípro- ca, temos o seguinte resultado. 46 Algoritmos Teorema 4.5. Seja N = (X,AX) uma network com n nós. A ultramétrica recíproca uRX pode ser calculada por uRX = ( max(AX , A T X) )(n−1) e a ultramétrica não-recíproca uNRX pode ser calculada por uNRX = max ( (AX)(n−1), (ATX)(n−1) ) , onde a operação (−)(n−1) é a (n− 1)-ésima potência da matriz na álgebra dioide. Demonstração. Por (4.4), temos A∗X = ũ∗X . Apenas trocando a notação, pela de�nição da ultramétrica não-recíproca temos que uNRX = max(A∗X , (A ∗ X)T ). Analogamente, a quasi-inversa da matriz simetrizada ĀX := max(AX , A T X) é [Ā∗X ]ij = min C(xi,xj) max k|xk∈C(xi,xj) ĀX(xk, xk+1). Pela de�nição de ultramétrica recíproca, segue que uRX = Ā∗X = max(AX , A T X)∗ e portanto, para concluir a demonstração, basta mostrar que A∗X = A (n−1) X . Note que na álgebra dioide U a operação ⊕ é idempotente, ou seja, a ⊕ a = min(a, a) = a, para todo a. Nesse caso, pela Proposição 4.3, temos que I ⊕ AX ⊕ A(2) X ⊕ · · · ⊕ A (k) X = (I ⊕ AX)(k), (4.5) para todo k ≥ 1. Além disso, como os elementos da diagonal são nulos nas matrizes do lado direito de (4.5) e os elementos fora da diagonal de I são +∞, segue que I ⊕ AX = AX . Consequentemente, I ⊕ AX ⊕ A(2) X ⊕ · · · ⊕ A (k) X = A (k) X . (4.6) Fazendo k → ∞ em ambos os lados da igualdade (4.6), e pela de�nição de quasi- inversa, obtemos A∗X = lim k→∞ A (k) X . Finalmente, pelo Teorema 4.4, segue que A(n−1) X = A (n) X , provando que o limite em (4.1) existe e, mais importante, que A∗X = A (n−1) X , como desejado. 5 Aplicações Neste capítulo, aplicaremos os métodos de clustering hierárquicos desenvolvidos anteriormente, ou seja, os métodos recíproco e não-recíproco, para obter dendrogramas para dois conjuntos de dados de networks assimétricas. Analisaremos a network de migração interna entre condados dos Estados Unidos (EUA) para o ano de 2011 e também analisaremos a network de migração entre países do mundo todo, no período 2000�2011. Para calcularmos computacionalmente as ultramétricas para esses dados, apresen- tamos inicialmente algumas de�nições. Denotamos por X o conjunto de condados dos Estados Unidos ou de países do mundo, e de�nimos a função µ : X × X → R+ ∪ {+∞}, em que µ(x, x′) é o número de indivíduos que migraram do condado (ou país) x para o condado (ou país) x′, e µ(x, x) = +∞, para todo x, x′ ∈ X. Assim, construímos uma network assimétrica NX = (X,AX), onde AX(x, x′) = { 0, se x = x′, f ( µ(x, x′)/ ∑ t µ(t, x′) ) , se x 6= x′, (5.1) para todo x, x′ ∈ X, onde f : [0, 1) → R+ é qualquer função decrescente. No nosso caso, optamos por considerar f(x) = 1− x. Analogamente, a função AX também pode ser de�nida tomando o somatório sobre a segunda variável, ou seja, AX(x, x′) = { 0, se x = x′, f ( µ(x, x′)/ ∑ t µ(x, t) ) , se x 6= x′. (5.2) A normalização µ(x, x′)/ ∑ t µ(t, x′) em (5.1) pode ser interpretada como sendo a probabilidade de que um imigrante do local x′ tenha vindo do local x, assim como µ(x, x′)/ ∑ t µ(x, t) em (5.2) pode ser interpretada como sendo a probabilidade de que um imigrante do local x tenha vindo do local x′. O papel da função f é transformar µ(x, x′)/ ∑ t µ(t, x′) e µ(x, x′)/ ∑ t µ(x, t) em dissimilaridades. 5.1 Migração da população dos EUA Iniciamos essa seção apresentando um resumo dos dados iniciais, a forma como estão e foram organizados. Tais dados estão no arquivo Net_Gross_US.txt obtido 47 48 Aplicações em [2]. A partir desse arquivo, �ltramos apenas as colunas 2, 3, 4, 5, 6, 8, produzindo um novo arquivo, raw_data_usa.csv, com 6 colunas, a saber: eA cA eB cB µ(cA, cB) µ(cB, cA) (5.3) onde e denota o nome de um dos estados1, c denota o nome de um condado do estado e, e µ é a migração (orientada) entre dois condados. O arquivo prepare_data_states_num_counties.csv contém informações sobre cada um dos 52 estados, mais precisamente, o nome do estado, sua quantia de condados, e seu código ISO, conforme a Tabela 5.1. O número total de condados é 3220, con- tando District of Columbia e Puerto Rico, e foram enumerados por ordem alfabética, primeiramente para os estados e então para os condados deste estado. Por exemplo, o primeiro condado do Alaska (segundo estado) recebe o número 68, pois existem 67 condados no único estado anterior, Alabama. Estado (em inglês) Cond. ISO Estado (em inglês) Cond. ISO Alabama 67 AL Montana 56 MT Alaska 29 AK Nebraska 93 NE Arizona 15 AZ Nevada 17 NV Arkansas 75 AR New Hampshire 10 NH California 58 CA New Jersey 21 NJ Colorado 64 CO New Mexico 33 NM Connecticut 8 CT New York 62 NY Delaware 3 DE North Carolina 100 NC District of Columbia 1 DC North Dakota 53 ND Florida 67 FL Ohio 88 OH Georgia 159 GA Oklahoma 77 OK Hawaii 5 HI Oregon 36 OR Idaho 44 ID Pennsylvania 67 PA Illinois 102 IL Puerto Rico 78 PR Indiana 92 IN Rhode Island 5 RI Iowa 99 IA South Carolina 46 SC Kansas 105 KS South Dakota 66 SD Kentucky 120 KY Tennessee 95 TN Louisiana 64 LA Texas 254 TX Maine 16 ME Utah 29 UT Maryland 24 MD Vermont 14 VT Massachusetts 14 MA Virginia 133 VA Michigan 83 MI Washington 39 WA Minnesota 87 MN West Virginia 55 WV Mississippi 82 MS Wisconsin 72 WI Missouri 115 MO Wyoming 23 WY Tabela 5.1: Conteúdo de prepare_data_states_num_counties.csv A seguir, explicaremos o conteúdo dos programas em Python criados para �ltrar os dados, extrair uma matriz de migração em condados, calcular as ultramétricas pelos 1Em todo o texto, também chamaremos de estado District of Columbia e Puerto Rico. Migração da população dos EUA 49 métodos recíproco e não-recíproco e, por �m, criar dendrogramas para cada network assim construída. Os programas em Python e os arquivos de dados estão disponíveis em [3]. 5.1.1 Preparando os dados: prepare_data.py Descreveremos abaixo os blocos mais relevantes do programa prepare_data.py. Vale ressaltar que as funções abaixo descritas estão em ordem alfabética, não na ordem em que serão executadas no programa. O programa inicia carregando os módulos necessários: NumPy e pandas, versões 1.13.3 e 0.19.2, respectivamente ([11, 9]). 1 # coding:utf-8 2 import numpy as np 3 import pandas as pd A função filter_line() recebe uma 6-upla como em (5.3), determina as posições i, j de cA, cB e de�ne o valor das entradas correspondentes na matriz de migração como sendo µ(cA, cB) e µ(cB, cA), ou seja, Mi,j = µ(cA, cB), Mj,i = µ(cB, cA). 4 def filter_line(s1, c1, s2, c2, x, y): 5 print 'filtering data to %s, %s, %s, %s' % (s1, c1, s2, c2) 6 from_ = (s1.decode('latin1'), c1.decode('latin1')) 7 to_ = (s2.decode('latin1'), c2.decode('latin1')) 8 from_idx = list_all_states_counties.index(from_) # defined in line 53 9 to_idx = list_all_states_counties.index(to_) # defined in line 53 10 migration_matrix[from_idx, to_idx] = y 11 migration_matrix[to_idx, from_idx] = x A função get_len() retorna o número de condados de um dado estado. 12 def get_len(state): 13 return int(len_states[np.where(len_states == state)[0][0], 1]) A função migration_from_states() é utilizada para obter os dados de migração a partir de um determinado estado, ou seja, �xado e ∈ states, a função cria dois arquivos: � mat_〈e〉_to_USA.csv: migração a partir de e para qualquer outro estado, não necessariamente diferente de e. Assim, representa uma matriz de ordemm×3220, onde m é o número de condados de e (obtido por get_len(〈e〉)) � mat_〈e〉_to_〈e〉.csv: migração interna apenas entre os condados de e, represen- tando uma matriz quadrada de ordem m, onde m é como acima. 14 def migration_from_states(): 15 migration_matrix = np.loadtxt(migration_matrix_file, dtype=int, delimiter=',')↪→ 16 for i, state in enumerate(states): 50 Aplicações 17 print 'filtering migration data from %s to USA' % state 18 len_state = get_len(state) 19 len_prev = sum([get_len(s) for s in states[:i]]) 20 tmp = migration_matrix[len_prev:len_prev+len_state] 21 tmp_single = migration_matrix[len_prev:len_prev+len_state, len_prev:len_prev+len_state]↪→ 22 np.savetxt('./migration_matrices/mat_%s_to_USA.csv' % state, tmp, fmt='%d', delimiter=',')↪→ 23 np.savetxt('./migration_matrices/mat_%s_to_%s.csv' % (state, state), tmp_single, fmt='%d', delimiter=',')↪→ A função migration_usa() recursivamente aplica a função filter_line() em cada linha do arquivo raw_data_usa.csv, salvando o resultado em migration_matrix.csv. Uma vez que os valores contidos nesse arquivo são os dados de migração entre dois condados, ele pode ser interpretado como uma matriz quadrada de ordem 3220. 24 def migration_usa(): 25 global migration_matrix 26 migration_matrix = np.zeros((num_counties, num_counties)) 27 with open(raw_data_usa_file, 'r') as f: 28 for line in f: 29 tmp = line.split(',') 30 if tmp[1] != '-' and tmp[3] != '-': 31 filter_line(tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5])↪→ 32 f.close() 33 np.savetxt(migration_matrix_file, migration_matrix, fmt='%d', delimiter=',')↪→ A função print_fields() é utilizada para �ltrar cada linha do Net_Gross_US.txt e obter os valores desejados, como em (5.3). 34 def print_fields(line): 35 data = '' 36 for cols in fields[used_fields]: 37 val = line[cols[0]-1:cols[1]].strip() 38 if val != '.': 39 data += val+',' 40 else: 41 data += val.replace('.', '')+',' 42 return data[:-1] A função use_pandas() utiliza o módulo Python pandas para ler o arquivo de dados raw_data_usa_file.csv, eliminar entradas não desejadas, como por exemplo, dados de migração de algum estado para outro país, ou de algum outro continente para algum estado, entre outros. Como resultado, temos a variável states que contém os 52 estados de interesse, bem como list_all_states_counties, que contém pares de estados e seus respectivos condados, ou seja, list_all_states_counties = {(e, c), c ∈ e, e ∈ states}. Migração da população dos EUA 51 Por �m, o número total de condados é calculado, a saber, 3220 e salvo na variável num_counties. 43 def use_pandas(): 44 global states, list_all_states_counties, num_counties 45 print 'Defining variables states, list_all_states_counties, num_counties'↪→ 46 df = pd.read_csv(raw_data_usa_file, encoding='latin1', header=None, delim_whitespace=False, index_col=None, sep=',', names=used_fields) ↪→ ↪→ 47 data = df[df[7] != '-'][used_fields] 48 states = data[4].unique().tolist() 49 states = states[:-9] 50 list_all_states_counties = [] 51 for state in states: 52 tmp = data[data[4] == state][5].unique().tolist() 53 list_all_states_counties += [(state, cty) for i, cty in enumerate(tmp)] #↪→ 54 num_counties = len(list_all_states_counties) A função write_raw_data() lê cada linha l do arquivo Net_Gross_US.txt e exe- cuta a função print_fields(l), escrevendo o resultado em uma linha do arquivo raw_data_usa.csv (que é utilizado pela função use_pandas()). 55 def write_raw_data(): 56 print 'Generating the file %s.\nWait...' % raw_data_usa_file 57 with open(raw_data_usa_file, 'w') as f: 58 with open(net_gross_file, 'rb') as net: 59 for l in net: 60 f.write('%s\n' % print_fields(l)) 61 net.close() 62 f.close() A seguir, apresentamos a parte principal do programa, na qual algumas variáveis são de�nidas e as funções acima são executadas. Primeiro, de�nimos as variáveis com os nomes dos arquivos que serão lidos (required �les) e que serão criados (output �les) por algumas funções. A variável used_fields informa quais colunas2 do arquivo Net_Gross_US.txt serão utilizadas. 63 """ MAIN CODE """ 64 65 """ Required files """ 66 net_gross_file = 'Net_Gross_US.txt' 67 fields = np.loadtxt('prepare_data_fields.csv', dtype=int, delimiter=',') 68 len_states = np.loadtxt('prepare_data_states_num_counties.csv', dtype=str, delimiter=',')↪→ 69 70 """ Output files """ 71 raw_data_usa_file = 'raw_data_usa.csv' 72 migration_matrix_file = 'migration_matrix.csv' 73 used_fields = [4, 5, 6, 7, 8, 10] 2A numeração das colunas inicia-se em 0. 52 Aplicações Por �m, executamos as funções abaixo para obtermos as matrizes de migração, que são salvas nos arquivos ./migration_matrices/*.csv. 74 write_raw_data() 75 use_pandas() 76 migration_usa() 77 migration_from_states() 78 79 """ End of file prepare_data.py """ 5.1.2 Criando as ultramétricas: compute_ultrametric.py O programa compute_ultrametric.py, como o nome diz, calcula as ultramétricas pelos métodos recíproco e não-recíproco, a partir das matrizes de migração, criadas no programa prepare_data.py da Seção 5.1.1. Lembramos que o algoritmo para calcular uR e uNR está descrito no Capítulo 3. Iniciamos o programa carregando o módulo NumPy. 1 # coding:utf-8 2 import numpy as np A função loop_ultrametric_inner_state() itera sobre todos os estados, exceto `District of Columbia', aplicando a função ultrametric_inner_state() em cada um deles que, de fato, calcula a ultramétrica para um estado dado. 3 def loop_ultrametric_inner_state(): 4 for state in states: 5 if state != 'District of Columbia': 6 ultrametric_inner_state(state) # defined in line 40 A função max_prod() implementa a operação de�nida em (4.3), que recebe duas ma- trizes quadradas de ordem n×n e calcula as entradas (i, j) e (j, i) da matriz resultante, por meio da função min_max(). Vale observar que, para o cálculo da ultramétrica, as duas matrizes de entrada mat1, mat2 são iguais. 7 def max_prod(mat1, mat2): 8 n = len(mat1) 9 tmp = np.zeros((n, n), dtype=float) 10 for i in range(n): 11 for j in range(i, n): 12 tmp[i, j] = min_max(mat1, mat2, i, j) 13 tmp[j, i] = min_max(mat2, mat1, j, i) 14 return tmp A função min_max() implementa a operação de�nida em (4.3). Graças ao módulo Numpy, essa operação pode ser feita simplesmente com uma linha de comando, usando min() e maximum(). 15 def min_max(mat1, mat2, i, j): 16 return min(np.maximum(mat1[i, :], mat2[:, j])) Migração da população dos EUA 53 A função power() recebe uma matriz A quadrada de ordem n e retorna a potência An−1, onde o produto é de�nido em (4.3). Devido à estabilidade do produto, como visto no Teorema 4.4, basta calcularmos as potências A2k até 2k ≤ n− 1. 17 def power(mat): 18 n = len(mat) 19 k = 1 20 while 2 ** (k - 1) < n: 21 mat = max_prod(mat, mat) 22 k += 1 23 return mat As duas seguintes funções, u_nr() e u_r() calculam as ultramétricas, pelos mé- todos não-recíproco e recíproco, respectivamente, a partir de uma matriz quadrada. Ambas funções possuem o argumento opcional output_file, cujos valores padrão são as strings 'uNR.csv' e 'uR.csv', que são os nomes dos arquivos .csv que conterão a matriz da ultramétrica e serão criados na pasta ./ultra_metrics. As funções implementam os cálculos descritos no Teorema 4.5. 24 def u_nr(mat, output_file='uNR.csv'): 25 print 'Computing non-reciprocal ultra metric. Wait... ', 26 tmp1 = power(mat) 27 tmp2 = tmp1.T 28 tmp = np.maximum(tmp1, tmp2) 29 np.savetxt('ultra_metric/'+output_file.replace(' ', '_'), tmp, fmt='%f', delimiter=',')↪→ 30 print 'saved in %s' % output_file.replace(' ', '_') 31 return tmp 32 33 def u_r(mat, output_file='uR.csv'): 34 print 'Computing reciprocal ultra metric. Wait... ', 35 a_bar = np.maximum(mat, mat.T) 36 tmp = power(a_bar) 37 np.savetxt('ultra_metric/'+output_file.replace(' ', '_'), tmp, fmt='%f', delimiter=',')↪→ 38 print 'saved in %s' % output_file.replace(' ', '_') 39 return tmp A função ultrametric_inner_state(), como o nome diz, trata da migração in- terna de um estado. Mais precisamente, dado um estado, carrega seu arquivo de migração interna e executa as funções u_r() e u_nr() nas duas matrizes retornadas pela função weight_fc(). Como consequência, as matrizes da ultramétrica são salvas nos arquivos .csv correspondentes. 40 def ultrametric_inner_state(state): # 41 print '\n** %s' % state 42 mat = np.loadtxt('migration_matrices/mat_%s_to_%s.csv' % (state, state), dtype=float, delimiter=',')↪→ 43 w_mat1, w_mat2 = weight_fc(mat) # defined in line 56 44 u_r(w_mat1, output_file='col_uR_%s.csv' % state) 45 u_nr(w_mat1, output_file='col_uNR_%s.csv' % state) 46 u_r(w_mat2, output_file='row_uR_%s.csv' % state) 54 Aplicações 47 u_nr(w_mat2, output_file='row_uNR_%s.csv' % state) A função ultrametric_national(), assim como a anterior, executa quatro vezes as funções que calculam as ultramétricas, porém, neste caso, a matriz é para o �uxo migratório nacional, ou seja, engloba também os dados entre diferentes estados. Vale observar que, devido à grande dimensão (32202) das matrizes, esta função pode levar horas para ser concluída. 48 def ultrametric_national(): 49 print 'This function needs time to finish.' 50 mat = np.loadtxt(migration_matrix_file, dtype=float, delimiter=',') 51 w_mat1, w_mat2 = weight_fc(mat) 52 u_r(w_mat1, output_file='col_uR.csv') 53 u_nr(w_mat1, output_file='col_uNR.csv') 54 u_r(w_mat2, output_file='row_uR.csv') 55 u_nr(w_mat2, output_file='row_uNR.csv') A função weight_fc() implementa os cálculos descritos em (5.1) e (5.2), para aplicar uma função peso à matriz de entrada. No código abaixo, duas funções são utilizadas, de modo que o processo retorna duas matrizes. As funções peso são f1(aij) = 1− aij/ ∑ k aik e f2(aij) = 1− aij/ ∑ k akj. 56 def weight_fc(mat): # 57 n = len(mat) 58 tmp1 = np.zeros((n, n), dtype=float) 59 tmp2 = np.zeros((n, n), dtype=float) 60 for j in range(n): 61 sum_col = max(mat[:, j].sum(), 1) 62 sum_row = max(mat[j, :].sum(), 1) 63 for i in range(n): 64 tmp1[i, j] = 1 - float(mat[i, j]) / sum_col 65 tmp2[i, j] = 1 - float(mat[i, j]) / sum_row 66 tmp1[j, j] = 0 67 tmp2[j, j] = 0 68 return tmp1, tmp2 Por �m, executamos as duas principais funções do programa, a partir do arquivo de dados de migração 'migration_matrix.csv'. 69 migration_matrix_file = 'migration_matrix.csv' 70 states = np.loadtxt('prepare_data_states_num_counties.csv', dtype=str, delimiter=',')[:, 0]↪→ 71 72 loop_ultrametric_inner_state() 73 74 """ This function could take hours to finish. """ 75 ultrametric_national() 76 77 """ End of file compute_ultrametric.py """ Migração da população dos EUA 55 5.1.3 Criando os dendrogramas: plot_dendrogram.py Nesta seção, apresentamos o último programa do projeto de migração nacional. O arquivo plot_dendrogram.py faz uso das matrizes de ultramétricas do diretório ultra_metric e cria no diretório dendrograms quatro dendrogramas para cada estado, e mais quatro para a migração nacional. Como de costume, carregamos os módulos necessários, dos quais Matplotlib [6] e Scipy [8] não haviam sido usados até então. O primeiro, fornece ferramentas para criarmos/exportarmos �guras e/ou grá�cos; o segundo, fornece funções para clustering. 1 # coding:utf-8 2 3 import matplotlib.pyplot as plt 4 import numpy as np 5 import pandas as pd 6 import sys 7 from scipy.cluster.hierarchy import dendrogram, linkage 8 9 sys.setrecursionlimit(2000) # to produce national dendrogram A função get_cmap() depende de n e para cada 0 ≤ i ≤ n retorna uma cor RGB distinta do mapa de cores name, cujo padrão é `hsv'. Tais cores serão usadas nas legendas dos dendrogramas. 10 def get_cmap(n, name='hsv'): 11 return plt.cm.get_cmap(name, n) A função get_iso() retorna o código ISO [7] de um dado estado, que será usado nas legendas dos dendrogramas. 12 def get_iso(state): 13 return len_states[np.where(len_states == state)[0][0], 2] A função get_len() retorna o número de condados em um dado estado. 14 def get_len(state): 15 return int(len_states[np.where(len_states == state)[0][0], 1]) A função loop_plot_dendrogram_by_state() itera sobre todos os estados (exceto `District of Columbia') e executa a função plot_dendrogram_by_state() neste estado, que gera os quatro dendrogramas para diferentes métodos e funções peso. 16 def loop_plot_dendrogram_by_state(): 17 for state in states: 18 if state != 'District of Columbia': 19 plot_dendrogram_by_state(state) # defined in line 51 A função plot_dendrogram_by_mat(), a partir de uma matriz de distâncias ou ultramétrica, aplica técnicas de clustering para criar um dendrograma. Os argumentos weight_fc e method são responsáveis pelas informações necessárias para criarmos o dendrograma, e seus possíveis valores são 'col' ou row', e 'R' ou 'NR', respectiva- mente. O argumento state recebe o nome do estado, se for o caso. 56 Aplicações 20 def plot_dendrogram_by_mat(mat, weight_fc='', method='', state='', label=np.array([])):↪→ 21 t = np.triu(mat, k=1) 22 tu = t[t > 0] 23 z = linkage(tu) 24 fig_width = 15 25 26 if len(mat) > 200: 27 fig_width = 25 28 else: 29 if len(mat) > 100: 30 fig_width = 20 31 if state == 'USA': 32 fig_width = 400 33 34 fig = plt.figure(figsize=(fig_width, 8)) 35 fig.add_subplot(111) 36 plt.title('Hierarchical Clustering Dendrogram of %s (wf: %s, meth: %s)' % (state, weight_fc, method))↪→ 37 dendrogram( 38 z, 39 color_threshold=0.9 * max(z[:, 2]), 40 leaf_rotation=90, 41 leaf_font_size=7, 42 labels=label, 43 ) 44 ax = plt.gca() 45 xlbls = ax.get_xmajorticklabels() 46 label_colors = {i: j for i, j in zip(label, colors)} 47 for lbl in xlbls: 48 lbl.set_color(label_colors[lbl.get_text()]) 49 plt.savefig('dendrograms/%s_u%s_%s.pdf' % (weight_fc, method, state), bbox_inches="tight")↪→ 50 plt.close(fig) A função plot_dendrogram_by_state() recebe um estado e itera sobre os mé- todos recíproco e não-recíproco e sobre as funções peso descritas em (5.1) e (5.2). Em cada iteração, carrega a matriz ultramétrica correspondente e executa a função plot_dendrogram_by_mat(). 51 def plot_dendrogram_by_state(state): # 52 print 'Creating all dendrograms for %s' % state 53 for method in ['R', 'NR']: 54 for weight_fc in ['col', 'row']: 55 mat = np.loadtxt('ultra_metric/%s_u%s_%s.csv' % (weight_fc, method, state.replace(' ', '_')), dtype=float, delimiter=',') ↪→ ↪→ 56 plot_dendrogram_by_mat(mat, weight_fc=weight_fc, method=method, state=state, label=dict_states_cts[state])↪→ Analogamente, a função plot_dendrogram_national() itera sobre os métodos re- cíproco e não-recíproco e sobre as funções peso e executa plot_dendrogram_by_mat() Migração da população dos EUA 57 para os dados de migração nacional. Vale observar que, devido à quantia de condados, a �gura do dendrograma precisa ser extremamente larga. Também, limitações de memória do computador podem interromper o procedimento. 57 def plot_dendrogram_national(): 58 print 'Creating the national migration dendrogram' 59 for method in ['R', 'NR']: 60 for weight_fc in ['col', 'row']: 61 mat = np.loadtxt('ultra_metric/%s_u%s.csv' % (weight_fc, method), dtype=float, delimiter=',')↪→ 62 plot_dendrogram_by_mat(mat, weight_fc=weight_fc, method=method, state='USA', label=all_counties)↪→ A função use_pandas(), como o nome sugere, utiliza o módulo pandas para gerar as variáveis dict_states_cts e states. A primeira é um dicionário Python, cuja chave é o nome de um estado e seu valor é uma lista com os condados deste estado. A segunda, é uma lista Python com os nomes dos estados. 63 def use_pandas(): 64 global dict_states_cts, states 65 print 'Defining variables dict_states_cts, states' 66 df = pd.read_csv(raw_data_usa_file, encoding='latin1', header=None, delim_whitespace=False, index_col=None, sep=',', names=used_fields) ↪→ ↪→ 67 data = df[df[7] != '-'][used_fields] 68 states = data[4].unique().tolist()[:-9] 69 70 dict_states_cts = {} 71 for state in states: 72 tmp = data[data[4] == state][5].unique() 73 tmp = [cty.replace(' County', '').replace(' Municipio', '').replace(' city', '').replace(' Parish', '') for cty in tmp] ↪→ ↪→ 74 dict_states_cts[state] = tmp Finalmente, após de�nidas as funções acima, o código a seguir carrega dois arquivos .csv, cria algumas variáveis (destacamos as listas Python all_counties e color) e executa as principais funções use_pandas(), loop_plot_dendrogram_by_state() e plot_dendrogram_national(). 75 raw_data_usa_file = 'raw_data_usa.csv' 76 len_states = np.loadtxt('prepare_data_states_num_counties.csv', dtype=str, delimiter=',')↪→ 77 used_fields = [4, 5, 6, 7, 8, 10] 78 79 use_pandas() 80 81 """ List of ISO codes + county names, 82 with commom words deleted to make them shorter """ 83 all_counties = [] 84 for state in states: 85 for county in dict_states_cts[state]: 58 Aplicações 86 tmp = u'%s %s' % (get_iso(state), county) 87 tmp = tmp.replace(' County', '').replace(' Municipio', '').replace(' city', '').replace(' Parish', '').replace(' Census Area', '') ↪→ ↪→ 88 all_counties.append(tmp) 89 90 colors = [] 91 cmap = get_cmap(len(states)) 92 for i, state in enumerate(states): 93 colors += [cmap(i)] * get_len(state) 94 95 # required files: 'ultra_metric/{col,row}_u{R,NR}_{state}.csv' 96 loop_plot_dendrogram_by_state() 97 98 # required files: 'ultra_metric/{col,row}_u{R,NR}.csv' 99 plot_dendrogram_national() 5.2 Migração da população mundial Nesta seção, faremos uma análise dos dados de migração mundial, criando ultramé- tricas e dendrogramas, como na seção anterior. O o programa world-migration.py é baseado nos descritos acima, com apenas alguns ajustes, pois os dados iniciais encontram-se em formato diferente. Assim, explicaremos somente os blocos de código que não foram já descritos. 5.2.1 Resumo dos dados iniciais A análise será feita a partir dos dados do arquivo worldbankdata.csv obtido em [4]. Este arquivo possui sete colunas, das quais �ltramos 1, 2, 5, 6, 7, contendo: país A cód. ISO país B cód. ISO �uxo de A para B O número de países é 231, separados em 7 regiões (que chamaremos de continentes, por abuso de terminologia), coloridas na legenda do dendrograma, do seguinte modo: AF Africa AS Asia CA Central America EU Europe NA North America OC Oceania SA South America Os dados de �uxo de A para B (e de B para A) são obtidos e guardados no arquivo migration_matrix.csv, que pode ser tratado como uma matriz quadrada de ordem 231. A seguir, com as mesmas funções do programa da seção anterior, as ultramétricas são calculadas para os dois métodos, recíproco e não-recíproco, e para as duas funções peso, como em (5.1) e (5.2). Por �m, quatro dendrogramas são criados. 5.2.2 Programa world-migration.py A seguir, apresentamos na íntegra o programa world-migration.py, apenas fa- zendo observações nos pontos relevantes e diferentes. Migração da população mundial 59 1 # coding:utf-8 2 3 import matplotlib.pyplot as plt 4 import numpy as np 5 import pandas as pd 6 from scipy.cluster.hierarchy import dendrogram, linkage A função delete_zero_col_row() é utilizada para eliminar linhas e/ou colunas nulas de uma matriz. 7 def delete_zero_col_row(mat): 8 idx = [] 9 for i in range(len(mat)): 10 if max(mat[i]) == 0: 11 idx.append(i) 12 13 mat = np.delete(mat, idx, 0) 14 mat = np.delete(mat, idx, 1) 15 return mat O bloco a seguir apresenta as funções já explicadas na seção anterior e portanto não requer mais detalhes. 16 def max_prod(mat1, mat2): 17 n = len(mat1) 18 tmp = np.zeros((n, n), dtype=float) 19 for i in range(n): 20 for j in range(i, n): 21 tmp[i, j] = min_max(mat1, mat2, i, j) 22 tmp[j, i] = min_max(mat2, mat1, j, i) 23 return tmp 24 25 def min_max(mat1, mat2, i, j): 26 return min(np.maximum(mat1[i, :], mat2[:, j])) 27 28 def plot_dendrogram_by_mat(mat, weight_fc='', method='', label=np.array([])):↪→ 29 t = np.triu(mat, k=1) 30 tu = t[t > 0] 31 z = linkage(tu) 32 fig_width = 15 33 34 if len(mat) > 200: 35 fig_width = 25 36 else: 37 if len(mat) > 100: 38 fig_width = 20 39 40 fig = plt.figure(figsize=(fig_width, 8)) 41 fig.add_subplot(111) 42 plt.title('Hierarchical Clustering Dendrogram (wf: %s, meth: %s)' % (weight_fc, method))↪→ 60 Aplicações 43 dendrogram( 44 z, 45 color_threshold=0.9 * max(z[:, 2]), 46 leaf_rotation=90, 47 leaf_font_size=7, 48 labels=label, 49 ) 50 ax = plt.gca() 51 xlbls = ax.get_xmajorticklabels() 52 label_colors = {i: j for i, j in zip(label, colors)} 53 for lbl in xlbls: 54 lbl.set_color(label_colors[lbl.get_text()]) 55 plt.savefig('%s_u%s.pdf' % (weight_fc, method), bbox_inches="tight") 56 plt.close(fig) 57 58 def plot_dendrogram_international(): 59 print 'Creating the international migration dendrogram' 60 for method in ['R', 'NR']: 61 for weight_fc in ['col', 'row']: 62 mat = np.loadtxt('%s_u%s.csv' % (weight_fc, method), dtype=float, delimiter=',')↪→ 63 plot_dendrogram_by_mat(mat, weight_fc=weight_fc, method=method, label=iso_codes)↪→ 64 65 def power(mat): 66 tmp = mat 67 n = len(tmp) 68 k = 1 69 while 2 ** (k - 1) < n: 70 tmp = max_prod(tmp, tmp) 71 k += 1 72 return tmp 73 74 def u_nr(mat, output_file='uNR.csv'): 75 print 'Computing non-reciprocal ultra metric. Wait... ', 76 tmp1 = power(mat) 77 tmp2 = tmp1.T 78 tmp = np.maximum(tmp1, tmp2) 79 np.savetxt(output_file.replace(' ', '_'), tmp, fmt='%f', delimiter=',')↪→ 80 print 'saved in %s' % output_file.replace(' ', '_') 81 return tmp 82 83 def u_r(mat, output_file='uR.csv'): 84 print 'Computing reciprocal ultra metric. Wait... ', 85 a_bar = np.maximum(mat, mat.T) 86 tmp = power(a_bar) 87 np.savetxt(output_file.replace(' ', '_'), tmp, fmt='%f', delimiter=',')↪→ 88 print 'saved in %s' % output_file.replace(' ', '_') Migração da população mundial 61 89 return tmp 90 91 def ultrametric_international(mat): 92 w_mat1, w_mat2 = weight_fc(mat) 93 u_r(w_mat1, output_file='col_uR.csv') 94 u_nr(w_mat1, output_file='col_uNR.csv') 95 u_r(w_mat2, output_file='row_uR.csv') 96 u_nr(w_mat2, output_file='row_uNR.csv') 97 98 def weight_fc(mat): 99 n = len(mat) 100 tmp1 = np.zeros((n, n), dtype=float) 101 tmp2 = np.zeros((n, n), dtype=float) 102 for j in range(n): 103 sum_col = max(mat[:, j].sum(), 1) 104 sum_row = max(mat[j, :].sum(), 1) 105 for i in range(n): 106 tmp1[i, j] = 1 - float(mat[i, j]) / sum_col 107 tmp2[i, j] = 1 - float(mat[i, j]) / sum_row 108 tmp1[j, j] = 0 109 tmp2[j, j] = 0 110 return tmp1, tmp2 O bloco a seguir de�ne várias variáveis, as mais importantes sendo raw_data_file e migration_matrix. As demais partes são detalhes técnicos, porém de fácil compre- ensão. 111 raw_data_file = 'worldbankdata.csv' 112 df = pd.read_csv(raw_data_file, header=0, delim_whitespace=False, index_col=0, sep=',')↪→ 113 iso_codes = df['Country Origin Code'].unique() 114 iso_codes_continents = np.loadtxt('codes.csv', dtype=str, delimiter=',') 115 iso_continents = sorted(list(set(iso_codes_continents[:, 1]))) 116 num_countries = len(iso_codes) 117 num_continents = len(iso_continents) 118 migration_matrix = df['value'].values.astype(float).reshape(num_countries, num_countries) ↪→ ↪→ 119 migration_matrix = delete_zero_col_row(migration_matrix) 120 np.savetxt('migration_matrix.csv', migration_matrix, fmt='%d', delimiter=',')↪→ 121 122 colors = [] 123 palette = ['orange', 'red', 'black', 'green', 'cyan', 'magenta', 'blue'] 124 for i, country in enumerate(iso_codes): 125 cont = iso_codes_continents[i, 1] 126 j = iso_continents.index(cont) 127 colors.append(palette[j]) Por �m, executamos as duas funções principais, que por sua vez, fazem uso das demais funções acima. A função ultrametric_international() recebe a variável migration_matrix e cria as ultramétricas de �uxo migratório internacional. A fun- 62 Aplicações ção plot_dendrogram_international(), como o nome já diz, cria arquivos .pdf dos quatro dendrogramas. 128 ultrametric_international(migration_matrix) 129 plot_dendrogram_international() 5.3 Ranking 5.3.1 Programa ranking_world_migration.py Neste programa, ainda sobre migração mundial, faremos uma análise para determi- nar quais países foram mais frequentes no �uxo migratório. Mais precisamente, dado um país P , �ltramos cinco países B1, . . . , B5 para os quais o �uxo de migração µ(P,Bi) de P para Bi é máximo, ou seja, �ltramos os cinco países que mais receberam população vinda de P . Analogamente, dualizamos o processo para obtermos A1, . . . , A5 para os quais o �uxo de migração µ(Ai, P ) de Ai para P é máximo. Um resumo dos dados obtidos pode ser visto na Tabela C.1 da página 81. Por �m, a partir da matriz de distâncias, criamos uma �ltração Fδ para δ ∈ { 7 10 , 8 10 , 9 10 }, onde Fδ = {(A,B); µ(A,B) ≤ δ}. Assim, podemos analisar se os países mais próximos com relação à migração também são mais próximos no sentido geográ�co, ou seja, se possuem fronteira comum ou se pertencem a um mesmo continente. Como de costume, carregamos os módulos necessários. 1 # coding:utf-8 2 3 import numpy as np 4 import pandas as pd A função filtration() utiliza a matriz de distâncias para criar a �ltração Fδ citada acima. Também, retorna o número de pares de países na �ltração que pertencem a um mesmo continente. 5 def filtration(mat): 6 delta_partition = [.7, .8, .9] 7 for delta in delta_partition: 8 same_continent = 0 9 idx, idy = np.where(mat <= delta) 10 pairs = zip(idx, idy) 11 for p in pairs: 12 tmp = '' 13 if p[0] < p[1]: 14 if iso_codes[p[0], 1] == iso_codes[p[1], 1]: 15 same_continent += 1 16 tmp = '** ' 17 print '%sdist %s: %s - %s %s-%s' % (tmp, mat[p[0], p[1]], countries[p[0]], countries[p[1]], iso_codes[p[0], 1], iso_codes[p[1], 1]) ↪→ ↪→ Ranking 63 18 print '%d/%d = %.02f belongs to same continent\n' % (same_continent, len(pairs)/2, float(same_continent)/(len(pairs)/2)) ↪→ ↪→ A função loop_top_migration_by_country() possui um laço que executa a função top_world_migration_by_country() para cada país country na variável countries. 19 def loop_top_migration_by_country(): 20 for country in countries: 21 print 'Ranking migration flow for %s' % country 22 top_world_migration_by_country(country) A função top_list() é interessante e requer uma explicação mais detalhada. Ao executarmos a função top_world_migration_by_country(), duas variáveis top_in_list e top_out_list são atualizadas, ou seja, cada uma recebe os cinco paí- ses de maior �uxo migratório com o país country em questão, como exempli�cado nas Figuras 5.1a e 5.1b. Após o laço realizado para cada país, as duas variáveis acima guardam todos os países que apareceram dentre os cinco primeiros. É aqui que a função top_list() faz o seu papel retornando, tanto para imigração quanto para emigração, os três países que mais aparecem dentre os top �ve. Informalmente, podemos pensar que, após todas as rodadas, são premiados nas ca- tegorias `Imigração' e `Emigração' com Ouro, Prata e Bronze os que mais frequentaram o pódio durante a competição (com o pódio sendo formado por cinco competidores). 23 def top_list(): 24 unique, counts = np.unique(top_in_list, return_counts=True) 25 print 'top in:\n', sorted(zip(counts, unique))[-3:] 26 unique, counts = np.unique(top_out_list, return_counts=True) 27 print 'top out:\n', sorted(zip(counts, unique))[-3:] A função top_world_migration_by_country(), para cada país country, deter- mina os (cinco) países com maior �uxo de migração para country e a partir de country. Ainda, salva os resultados nos arquivos ranking/top_from_.csv e ranking/top_to_.csv. 28 def top_world_migration_by_country(country, save_csv=False): 29 from_ = data[data['from'] == country] 30 top_from = from_.sort_values(['value'], ascending=False).head(5)[['from', 'to', 'value']]↪→ 31 to_ = data[data['to'] == country] 32 top_to = to_.sort_values(['value'], ascending=False).head(5)[['to', 'from', 'value']]↪→ 33 34 for val in top_from['to'].values: 35 top_in_list.append(val) 36 for val in top_to['from'].values: 37 top_out_list.append(val) 38 39 if save_csv: 64 Aplicações 40 top_from.to_csv('ranking/top_from_%s.csv' % country, sep=',', encoding='utf-8', index=False, header=False)↪→ 41 top_to.to_csv('ranking/top_to_%s.csv' % country, sep=',', encoding='utf-8', index=False, header=False)↪→ A função top_world_migration() simplesmente retorna os top ten, ou seja, os dez países com maior �uxo de imigração e emigração. Ainda, salva os resultados nos arquivos ranking/top_out.csv e ranking/top_in.csv. 42 def top_world_migration(save_csv=False): 43 """ This function saves the top 10 countries 44 with higher incoming/outcoming flow. """ 45 top_out = data.groupby(['from'], as_index=False).sum().sort_values(['value'], ascending=False).head(10)[['from', 'value']].astype(object) ↪→ ↪→ 46 top_in = data.groupby(['to'], as_index=False).sum().sort_values(['value'], ascending=False).head(10)[['to', 'value']].astype(object) ↪→ ↪→ 47 if save_csv: 48 top_out.to_csv('ranking/top_out.csv', sep=',', encoding='utf-8', index=False, header=False)↪→ 49 top_in.to_csv('ranking/top_in.csv', sep=',', encoding='utf-8', index=False, header=False)↪→ A seguir, consta a parte principal do programa, onde o arquivo worldbankdata.csv de dados bruto ([4]) é carregado e as funções descritas acima são executadas. Tam- bém, é carregado o arquivo row_uNR.csv (ou os outros três possíveis) da matriz da ultramétrica, que é utilizada pela função filtration(). 50 """ Source: Global Bilateral Migration 51 Data from 2000 52 Last Updated: 06/28/2011 """ 53 raw_data_file = 'worldbankdata.csv' 54 col_names = ['from', 'iso_from', 'tmp1', 'tmp2', 'to', 'iso_to', 'value'] 55 df = pd.read_csv(raw_data_file, encoding='utf8', header=0, delim_whitespace=False, index_col=None, sep=',', names=col_names)↪→ 56 data = df[['from', 'iso_from', 'to', 'iso_to', 'value']] 57 countries = df['from'].unique().tolist() 58 iso_codes = np.loadtxt('codes.csv', delimiter=',', dtype=str) 59 top_in_list = [] 60 top_out_list = [] 61 62 top_world_migration() 63 loop_top_migration_by_country() 64 top_list() 65 66 mat = np.loadtxt('row_uNR.csv', delimiter=',') 67 filtration(mat) Ranking 65 5.3.2 Classi�cando os dados: ranking_USA_migration.py Neste programa, faremos a mesma análise da Seção 5.3.1 para determinar quais estados foram mais frequentes no �uxo migratório. Nesse caso, dado um estado S, �ltramos cinco estados E1, . . . , E5 para os quais o �uxo de migração µ(S,Ei) de S para Ei é máximo. Analogamente, dualizamos o processo para obtermos E1, . . . , E5 para os quais o �uxo de migração µ(Ei, S) de Ei para S é máximo. Como o nosso banco de dados raw_data_usa.csv ([2]) é muito maior, não faremos a �ltração para analisar se os estados próximos com relação à migração também são mais próximos no sentido geográ�co. Primeiro, carregamos os módulos necessários. 1 # coding:utf-8 2 3 import numpy as np 4 import pandas as pd A função loop_top_migration_by_state() é análoga à função da seção anterior loop_top_migration_by_country(). 5 def loop_top_migration_by_state(): 6 for state in states: 7 print 'Ranking migration flow for %s' % state 8 top_migration_by_state(state) A função top_list() retorna os três estados de maior �uxo migratório. 9 def top_list(): 10 unique, counts = np.unique(top_from_list, return_counts=True) 11 print 'top out:\n', sorted(zip(counts, unique))[-3:] 12 unique, counts = np.unique(top_to_list, return_counts=True) 13 print 'top in:\n', sorted(zip(counts, unique))[-3:] A função top_migration_by_state() é análoga à função top_world_migration_ by_country() e salva os resultados nos arquivos ranking/top_from_.csv e ranking/top_to_.csv. 14 def top_migration_by_state(state, save_csv=False): 15 tmp = data[data[4] == state] 16 top_from = tmp.groupby([6], as_index=False).sum().sort_values([10], ascending=False).head(5)[[6, 10]]↪→ 17 top_to = tmp.groupby([6], as_index=False).sum().sort_values([8], ascending=False).head(5)[[6, 8]]↪→ 18 19 for val in top_from[6].values: 20 top_from_list.append(val) 21 for val in top_to[6].values: 22 top_to_list.append(val) 23 24 if save_csv: 25 top_from.to_csv('ranking/top_from_%s.csv' % state, sep=',', encoding='utf-8', index=False, header=False)↪→ 66 Aplicações 26 top_to.to_csv('ranking/top_to_%s.csv' % state, sep=',', encoding='utf-8', index=False, header=False)↪→ S E1 E2 E3 E4 E5 (a) top_to_list S E1 E2 E3 E4 E5 (b) top_from_list Figura 5.1: Representação da função top_migration_by_state() A função top_migration_national() retorna os top ten dos estados com maior �uxo migratório, listados na Tabela C.3 da página 82. Ainda salva os resultados nos arquivos ranking/top_in_USA.csv e ranking/top_out_USA.csv. 27 def top_migration_national(save_csv=False): 28 print 'Ranking national migration...', 29 top_in = data.groupby([6, 7], as_index=False).sum().sort_values([10], ascending=False).head(10)[[6, 7, 10]]↪→ 30 top_out = data.groupby([6, 7], as_index=False).sum().sort_values([8], ascending=False).head(10)[[6, 7, 8]]↪→ 31 if save_csv: 32 top_in.to_csv('ranking/top_in_USA.csv', sep=',', encoding='utf-8', index=False, header=False)↪→ 33 top_out.to_csv('ranking/top_out_USA.csv', sep=',', encoding='utf-8', index=False, header=False)↪→ 34 print 'done.' Finalmente, a parte principal do programa, onde os arquivos raw_data_usa.csv e prepare_data_states_num_counties.csv são carregados e as funções acima são executadas. 35 """ MAIN CODE """ 36 37 raw_data_usa_file = 'raw_data_usa.csv' 38 states = np.loadtxt('prepare_data_states_num_counties.csv', dtype=str, delimiter=',')[:, 0]↪→ 39 40 used_fields = [4, 5, 6, 7, 8, 10] 41 df = pd.read_csv(raw_data_usa_file, encoding='latin1', header=None, delim_whitespace=False, index_col=None, sep=',', names=used_fields)↪→ 42 data = df[(df[5] != '-') & (df[7] != '-')][used_fields] 43 Ranking 67 44 top_from_list = [] 45 top_to_list = [] 46 47 top_migration_national() 48 loop_top_migration_by_state() 49 top_list() Referências [1] G. Carlsson, F. Mémoli, A. Ribeiro, and S. Segarra. Axiomatic construction of hierarchical clustering in asymmetric networks. Acoustics, Speech and Signal Pro- cessing (ICASSP), 2013 IEEE International Conference, pages 5219�5223, 2013. [2] County-to-county migration �ows: 2011-2015 acs, 2018. [3] T. de Melo and B. Barreiro. Python program to compute ultra-metric for migration networks, 2018. [4] Global Bilateral Migration Database, 2018. [Online; acessado em 14 de dezembro de 2018]. [5] M. Gondran and M. Minoux. Graphs, dioids and semi rings: New models and algorithms. Springer, New York, 2008. [6] J. D. Hunter. Matplotlib: A 2d graphics environment. Computing In Science & Engineering, 9(3):90�95, 2007. [7] Iso 3166-1 � Wikipedia, the free encyclopedia, 2018. [Online; acessado em 14 de dezembro de 2018]. [8] E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scienti�c tools for Python, 2001�. [Online; acessado em 14 de dezembro de 2018]. [9] W. McKinney. Data structures for statistical computing in python. Proceedings of the 9th Python in Science Conference, 2010�. [Online; acessado em 14 de dezembro de 2018]. [10] Migracão � Wikipedia, the free encyc