Embed Size (px)
Citation preview
Universidade de BrasíliaInstituto de Ciências Exatas
Departamento de Ciência da Computação
Algoritmos Multiobjetivos para PlanejamentoSistemático de Conservação: Estudo de Caso para
Plantas e Mamíferos Terrestres
Gustavo Fernandes de Almeida
Monografia apresentada como requisito parcialpara conclusão do Curso de Engenharia da Computação
OrientadoraProf.ª Dr.ª Maria Emília Machado Telles Walter
Brasília2016
Universidade de BrasíliaInstituto de Ciências Exatas
Departamento de Ciência da Computação
Algoritmos Multiobjetivos para PlanejamentoSistemático de Conservação: Estudo de Caso para
Plantas e Mamíferos Terrestres
Gustavo Fernandes de Almeida
Monografia apresentada como requisito parcialpara conclusão do Curso de Engenharia da Computação
Prof.ª Dr.ª Maria Emília Machado Telles Walter (Orientadora)CIC/UnB
Prof.ª Dr.ª Célia Ghedini Ralha Prof. Dr. Guilherme Novaes RamosCIC/UnB CIC/UnB
Prof. Dr. Ricardo Pezzuol JacobiCoordenador do Curso de Engenharia da Computação
Brasília, 1º de julho de 2016
Dedicatória
A minha família, mesmo sem saber ao certo o que eu faço, me apoiouincondicionalmente.
A todos aqueles que perguntaram como Julia está.
iii
Agradecimentos
A professora Maria Emília, pela orientação e paciência.A professora Célia Ghedini Ralha, pelos seus comentários precisos e ao professor Guil-
herme Novaes Ramos, pelas observações francas.A Drª. Shana Schlottfeldt Santos, Dr. Ricardo Dobrovolski e a Drª. Mariana Pires de
Campos Telles, pela disposição em esclarecer as diversas dúvidas que surgiram duranteesse trabalho.
Finalmente, aos amigos.
iv
Resumo
Este trabalho apresenta a implementação em linguagem Julia de método para indicaçãode áreas de conservação ambiental baseado em algoritmos de otimização multiobjetivo.O método dá suporte à decisão para o planejamento de novas áreas de conservação,considerando simultaneamente e sem preferências diversos objetivos antagônicos, quaissejam, minimização de áreas a serem preservadas, considerando diversidade genética eequilíbrio de Hardy-Weinberg. Verificou-se sua aplicabilidade em duas diferentes situ-ações, primeiramente replicando um experimento anterior feito sobre a planta do CerradoBrasileiro Baru para assegurar a correção da presente implementação e depois uma apli-cação considerando mamíferos terrestres no Distrito Federal e na América do Sul. Osresultados mostram que o método gera um portfólio de alternativas que consideram si-multaneamente objetivos conflitantes, que podem ser usados como ponto de partida pelosresponsáveis pelas tomadas de decisão para auxiliar suas escolhas.
Palavras-chave: Ecoinformática, Otimização Multiobjetivos, NSGA-II, linguagem Julia,Planejamento Sistemático de Conservação, Cerrado Brasileiro, Plantas, Baru, Mamíferos
v
Abstract
This project details a method’s implementation to help choose conservation areas basedon multi-objective optimization. This method aims to help decision makers providinga set of possible solutions that simultaneously ,and without preference, consider manyantagonistic objectives, those may be minimizing conservation areas, considering geneticdiversity and Hardy-Weinberg equilibrium. The applicability of the method was evaluatedin two situations, redoing a previous experiment to assure the implementation’s correct-ness and after considering terrestrial mammals on the Brazilian Federal District and onSouth America. The method’s results showed a portfolio of alternatives that consideredconflicting goals that could be used as starting point to choose new conservation areas.
Keywords: Ecoinformatics, Multiobjective optimization, NSGA-II, Julia language, Mam-mals, Plants, Baru, Brazilian Cerrado, Systematic Conservation Planning
vi
Sumário
1 Introdução 11.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Descrição dos Capítulos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Planejamento Sistemático de Conservação 52.1 Conceitos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Ferramentas computacionais . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Marxan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.2 C-Plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2.3 TARGET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Otimização Multiobjetivo para PSC . . . . . . . . . . . . . . . . . . . . . . 11
3 Otimização 143.1 Conceitos Básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1 Otimização Multiobjetivo . . . . . . . . . . . . . . . . . . . . . . . 143.1.2 Frente de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Algoritmos Evolutivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2.1 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2.2 Algoritmos Evolutivos e Otimização Multiobjetivo . . . . . . . . . . 19
3.3 Algoritmos Genéticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.4 Algoritmo NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.1 NSGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.4.2 Descrição do NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 NSGA-II usando Julia 334.1 Otimização usando Julia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.1 Linguagem Julia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.2 Detalhes da implementação do NSGA-II em Julia . . . . . . . . . . . . . . 35
vii
5 Estudo de Caso: Baru 385.1 Baru . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385.2 Dados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2.1 Matrizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.3 Simulações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.1 Simulação 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 495.3.2 Simulações 4D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.3.3 Simulação 5D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6 Estudo de Caso: Mamíferos Terrestres 576.1 Dados dos mamíferos terrestres . . . . . . . . . . . . . . . . . . . . . . . . 576.2 Presença e ausência de mamíferos terrestres . . . . . . . . . . . . . . . . . 57
6.2.1 Distrito Federal e Entorno . . . . . . . . . . . . . . . . . . . . . . . 576.2.2 América do Sul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7 Conclusões e Trabalhos Futuros 67
Referências 68
Anexo 71
I Código-fonte do NSGA-II em Julia 72
II Matriz A 92
IIIMatriz B 94
viii
Lista de Figuras
1.1 Mapa com as áreas prioritárias para a conservação no bioma Cerrado ePantanal1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Zonae Cogito com Marxan e GIS. . . . . . . . . . . . . . . . . . . . . . . . 82.2 Distribuição das populações estudadas [31]. . . . . . . . . . . . . . . . . . . 12
3.1 Representações gráficas para otimização em duas dimensões, na qual busca-se minimizar ambas as dimensões. . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Como os algoritmos evolutivos pode ser classificados. . . . . . . . . . . . . 183.3 Fluxograma ilustrando algoritmo evolutivo simples, adaptada do Capítulo 2
de Analyzing Evolutionary Algorithms [20]. As setas contínuas representamo fluxo de controle e as setas tracejadas o fluxo de dados. . . . . . . . . . . 20
3.4 Gráfico da Equação 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.5 Representação de um indivíduo na população. . . . . . . . . . . . . . . . . 233.6 Figura esquemática de um torneio binário. . . . . . . . . . . . . . . . . . . 243.7 Recombinação de ponto único. . . . . . . . . . . . . . . . . . . . . . . . . . 253.8 Mutação bit-flip no gene 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.9 Um conjunto de soluções para otimização 5D, resultante de uma execução.
Visão do número de populações selecionadas, número de alelos faltantes efrequência alélica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.10 As etapas do algoritmo NSGA-II, sendo as etapas indicadas de 1 a 6. . . . 30
4.1 Gráfico do benchmark comparando Julia com outras linguagens em escalalogarítmica [6]. O eixo das ordenadas indica é o tempo decorrido em com-paração à implementação em C e o das abscissas são os diversos algoritmosimplementados em várias linguagens. . . . . . . . . . . . . . . . . . . . . . 35
5.1 Árvore e fruto do Dipteryx alata. . . . . . . . . . . . . . . . . . . . . . . . 395.2 Representação do vetor solução, com índices relativos às e as respectivas
populações locais que representam. Neste exemplo, as populações escolhi-das foram CMT, PCO, AMS, ATO, MAMG, STGO, AMG, CMS e ARTO. 49
ix
5.3 Um conjunto de soluções para uma execução do NSGA-II para duas di-mensões. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.4 Um conjunto de soluções para uma execução do NSGA-II para duas di-mensões. Soluções que ficariam representadas de maneira sobreposta foramdeslocadas para melhor visualização. . . . . . . . . . . . . . . . . . . . . . 51
5.5 Um conjunto de soluções para otimização 5D, resultante de uma execução.Visão do número de populações selecionadas, número de alelos faltantes efrequência alélica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.6 Um conjunto de soluções para otimização 5D, resultante de uma execução.Visão do número de populações selecionadas, valor de heterozigose e EHW. 56
6.1 Área analisada durante o experimento no DF e Entorno. . . . . . . . . . . 606.2 As soluções encontradas para otimização 2D no DF e entorno. . . . . . . . 636.3 A América do Sul, com quatro pontos extremos em destaque, e área recortada. 646.4 Regiões a serem preservadas na América do Sul, com coordenadas de con-
servação destacadas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
x
Lista de Tabelas
5.1 Nome das populações amostradas e quantidade de indivíduos amostrados. 405.2 Alelos identificados para os 9 loci sequenciados. . . . . . . . . . . . . . . . 415.3 Representação do parcial dos dados utilizados . . . . . . . . . . . . . . . . 415.4 Valores de homozigose e heterozigose para cada população analisada. . . . 425.5 Matriz intermediária parcial X que relaciona indivíduos com os alelos iden-
tificados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.6 Representação parcial da Matriz A - relaciona populações com os alelos
identificados. A Matriz A transposta na íntegra pode ser encontrada noAnexo II. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.7 Representação parcial da Matriz B - frequência dos alelos dentro das po-pulações locais normalizada para quantidade de amostras por população.Valores multiplicados por 104 para melhor visualização. A Matriz B trans-posta na íntegra pode ser encontrada no Anexo III. . . . . . . . . . . . . 45
5.8 Representação parcial da primeira matriz intermediária para elaboração daMatriz C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.9 Representação parcial da segunda matriz intermediária para elaboração daMatriz C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.10 Matriz C - associa os valores de heterozigose dos alelos identificados compopulações locais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.11 Valores de EHW de cada população analisada. . . . . . . . . . . . . . . . . 485.12 Resultados encontrados para nenhum alelo faltante com as sete populações
selecionadas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 505.13 Resultados de irreplaceability para nenhum alelo faltante com as sete po-
pulações selecionadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.14 Tabela que agrega os resultados encontrados para duas otimizações 4D com
7 e 8 populações. Todos as soluções apresentadas englobam todos alelos. . 545.15 Resultados das soluções para otimização 5D para 7 e 8 populações. . . . . 55
xi
6.1 Representação parcial da matriz de presença/ausência de mamíferos ter-restres em determinadas regiões geográficas. As coordenadas geográficassão dadas em notação decimal. . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.2 Tabela parcial de valores de presença e ausência de mamíferos terrestres decada par de coordenadas do DF e Entorno. . . . . . . . . . . . . . . . . . 59
6.3 Conjunto das soluções encontradas para otimização 2D na área definidapara o DF e entorno com 5 e 6 populações selecionadas. Em todos osresultados, todos os mamíferos são incluídos. . . . . . . . . . . . . . . . . . 61
6.4 Tabela parcial da MAS com valores de presença e ausência de mamíferosterrestres de cada par de coordenadas da América do Sul. . . . . . . . . . 65
II.1 Matriz A transposta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
III.1 Matriz B transposta. Valores multiplicados por 104 para melhor visualização. 95
xii
Lista de Abreviaturas e Siglas
EHW Equilíbrio de Hardy-Weinberg.
GIS Geographic Information System.
JIT Just-in-time.
MOGA Multi-Objective Genetic Algorithm.
NSGA Nondominated Sorting Genetic Algoritm.
PAES Pareto Archived Evolution Strategy.
PSC Planejamento Sistemático de Conservação.
SPEA Strenght Pareto Evolutionary Algorithm.
VEGA Vector Evaluated Genetic Algorithm.
xiii
Capítulo 1
Introdução
Com o crescimento das cidades, sobretudo no Brasil, demandas por produtos e servi-ços aumentaram [8, 21]. Essa realidade impõe a busca por métodos de produção dematérias-primas capazes de atender essas necessidades. Em particular, a industrializa-ção da atividades agrícolas se impuseram, provocando uma forte expansão das fronteirasagrícolas no país.
No centro das discussões de como melhor utilizar os recursos naturais está a questão dapreservação ambiental. A necessidade estratégica [36] da preservação do capital natural,expressa na biodiversidade e nas diversas ações prestados pelo meio ambiente como apolinização parcialmente responsável pelo aumento da produção de plantas ou as matasque previnem deslizamento de terra, deve ser discutida com prioridade.
Maneiras de identificar e definir quais são as áreas cuja preservação é imprescindívelpara manter a biodiversidade tornam-se cada vez mais importantes. Esse fato entra emdestaque em países em desenvolvimento [7, 21], visto que nesses países a preservaçãonão é prioridade econômica, enquanto o ganho financeiro a curto prazo com a utilizaçãodesordenada das possíveis áreas de preservação é economicamente muito mais atrativa.
A escolha das áreas a serem preservadas não deve ser feita de maneira aleatória.Com incentivos à apropriação de áreas ainda não cultivadas pela indústria agrícola, ea busca de recursos para a criação de reservas, a escolha dessas novas reservas deveser cuidadosamente realizada para obter o melhor custo-benefício entre a expansão dasfronteiras agrícolas e a preservação das espécies nativas.
Infelizmente o processo de escolha de uma reserva natural não é simples. As questõestécnicas, relacionadas com os campos da Biologia e Ecologia, compõem apenas uma fraçãode todas os tópicos a serem considerados. Essas considerações práticas são essenciais paraa validação, relevância e a capacidade da reserva atingir seus objetivos preservacionistas,mas a criação de uma nova reserva apresenta viés político, econômico e social, que sãoespecialmente complicados de mensurar.
1
Os aspectos que fogem do escopo da Ecologia são tão importantes quanto os que elaengloba, no sentido que o projeto de uma reserva não é suficiente para a preservaçãoambiental. Questões socioeconômicas como a maneira de implementar e manter umareserva são indispensáveis e tão importantes quanto outros requisitos técnicos.
Um questionamento que surge é como harmonizar esse conjunto de demandas tãoheterogêneo, respeitando as exigências e restrições de cada uma? Nesse projeto, buscamosaprofundar a discussão de como obter respostas para essa pergunta, usando o métodoproposto de otimização multiobjetivo usando informações de plantas a nível genérico [31],para apontar locais adequados de preservação para uma espécie de planta do CerradoBrasileiro.
A Figura 1.1 mostra o mapa definido na Portaria Ministério do Meio Ambiente n.º09/20071, essa Portaria fez um levamento de áreas prioritárias para conservação, e bene-fícios da biodiversidade brasileira do Cerrado e Pantanal ilustra a necessidade de preser-vação desse bioma. O levantamento coloca 431 áreas como prioritárias, das quais apenas181 são áreas protegidas (unidades de conservação e terras indígenas).
Figura 1.1: Mapa com as áreas prioritárias para a conservação no bioma Cerrado ePantanal1.
1http://mma.gov.br/biomas/cerrado/areas-prioritarias, acesso 3 de julho de 2016.
2
1.1 Motivação
No cenário brasileiro atual, as iniciativas para conservação ambiental são carentes de meca-nismos que auxiliem a difícil tarefa de equilibrar diferentes agendas de grupos interessados.Essa falta é agravada pela limitada variedade de soluções desenvolvidas especificamentepara as necessidades da realidade vivenciada no Brasil.
Além disso, em nível mundial, dados de mamíferos terrestres foram obtidos e ini-cialmente analisados, usando algoritmos mono-objetivos e depois pós-processados paracontemplar os objetivos conflitantes.
1.2 Problema
O problema que buscamos solucionar é verificar se existem métodos computacionais paradeterminar um conjunto de soluções que atendam da melhor forma possível vários ob-jetivos conflitantes, para indicar reservas ambientais. Esses objetivos incluem aspectosespecíficos das espécies estudadas e outros interesses diversos como ocupação humana,agricultura e disponibilidade de água.
1.3 Objetivos
O objetivo geral desse projeto é utilizar algoritmo de otimização multiobjetivo para dedeterminar áreas de conservação ambiental considerando características específicas deplantas e animais mamíferos terrestres.
Os objetivos específicos são:
• Implementar o algoritmo NSGA-II, de otimização multiobjetivo;
• Avaliar o NSGA-II com duas, quatro e cinco dimensões, para indicar áreas de pre-servação da planta Dypteryx alata (baru). As dimensões representam presença/au-sência de alelos em uma população, número de populações, frequência de alelos,heterozigose e equilíbrio de Hardy-Weinberg (EHW);
• Usar o NSGA-II, em duas dimensões, para indicar áreas de preservação no DistritoFederal e entorno, e também na América do Sul, para mamíferos terrestres. Osdados indicam a presença e ausência de mamíferos no mundo todo, neste trabalho,utilizamos apenas as informações relativas a essas duas regiões.
3
1.4 Descrição dos Capítulos
O Capítulo 2 trata do Planejamento Sistemático de Conservação, uma área que defineetapas para o processo de conservação da natureza. Em seguida apresenta algumas ferra-mentas disponíveis na literatura para definição de áreas de conservação. Depois, trata docaso específico estudado por Santos et al. [31], com detalhes da abordagem e dos resultadosobtidos.
O Capítulo 3 traz conceitos importantes em otimização. Primeiramente apresentaconceitos de otimização multiobjetivo, algoritmos evolutivos e algoritmos genéticos. Aofinal, discorre sobre os algoritmos NSGA e NSGA-II.
O Capítulo 4 apresenta brevemente a linguagem Julia e mostra detalhes da implemen-tação do NSGA-II em Julia.
O Capítulo 5 discute primeira série de experimentos desse trabalho, com o objetivo dereplicar aquele desenvolvido em Santos et al. [31] para a planta do Cerrado Dypteryx alata,o baru. Os dados utilizados são descritos, assim como esses dados foram trabalhados e assimulações feitas.
O Capítulo 6 apresenta a segunda série de experimentos, que tratam da otimizaçãoda presença de mamíferos terrestres na área do Distrito Federal e entorno e em seguidaem toda a América do Sul.
O Capítulo 7 conclui o trabalho e apresenta proposta de trabalhos futuros.
4
Capítulo 2
Planejamento Sistemático deConservação
Neste capítulo, apresentamos conceitos básicos sobre a Planejamento Sistemático de Con-servação (PSC), necessários ao entendimento deste trabalho. Buscamos explicar o que éo processo Planejamento Sistemático de Conservação e suas etapas; apresentar ferramen-tas amplamente utilizadas para o auxílio da definição de áreas e descrever o estudo quenorteia esta abordagem.
2.1 Conceitos básicos
A necessidade de proteger territórios sempre foi foco de discussão territorial das nações.No seu núcleo, estão perspectivas e prioridades de como melhor aproveitar o espaço, sejamantendo-o intacto, modificando-o fortemente ou qualquer outra opção intermediária.Com o objetivo de nortear e melhor embasar debates sobre conservação, foram compiladasestratégias para efetivamente atingir metas, incluindo pontos de vista antagônicos. O PSC(Systematic Conservation Planning) define uma série de etapas para indicar as melhoresdecisões para conservação [26].
Reservas naturais surgem espontaneamente quando falamos de preservação ambien-tal. Dois são os principais papéis de reservas [26]: devem representar a biodiversidadede alguma região e devem separar essa região de qualquer atividade que ameace a suapermanência. Sem dúvida elas representam uma importante parte no esforço de conserva-ção, mas também são tendenciosas. Elas representam, em grande parte, locais intocáveise com baixas capacidades comerciais. Portanto não devem ser vistas como vias exclusivasde preservação. O PSC leva em consideração esse viés durante seu desenvolvimento.
O PSC é dividido em seis passos [26], que não são necessariamente sequenciais, ouseja, pode haver retroalimentação entre diferentes etapas, e são mutáveis, podem sofrer
5
modificações para adequação a diferentes contextos. O primeiro passo é o levantamentode dados sobre a biodiversidade da região planejada, seguido pela definição de objetivosclaros para a conservação da região. Depois temos a revisão das áreas de conservação jáestabelecidas, seleção de novas áreas, implementação de ações conservativas e, por fim, amanutenção dos valores de conservação determinados para elas.
A primeira etapa do PSC é de preparação e a ideia é fazer amostragem da biodiver-sidade, suficiente para representar com algum grau de confiança a região como um todo.Não é uma tarefa trivial e, normalmente, depende bastante de dados coletados. Caso setenham recursos suficientes, coletam-se mais dados para ampliar ou substituir os dadospresentes. Espécies raras ou ameaçadas na região devem receber atenção especial. Porexemplo, pode-se utilizar otimização multiobjetivo para ajudar no projeto de criação deuma coleção que preserve plantas ou animais, denominado germoplasma [32].
A segunda etapa é identificar objetivos de conservação para a região. O foco é trazermetas, preferencialmente quantitativas, para uso operacional. Para definir tais metasdiferentes estratégias são adotadas como a inclusão de stakeholders1, tanto especialistasquanto leigos para aumentar a influência política da conservação ambiental [29], o usode simulações de modelos climáticos [11], ou o estudo da densidade da população dedeterminadas espécies [39].
A terceira etapa visa fazer o levantamento do estado atual da região. Procuram-se asfalhas nas redes de reservas, que mostram os aspectos bem supridos e quais necessitam deatenção. Por outro lado, o uso exclusivo da análise por falhas pode causar complicaçõesa longo prazo, visto que esse tipo de análise é imediatista e tem dificuldade em analisar adinâmica e persistência da biodiversidade no futuro.
A quarta etapa engloba selecionar reservas adicionais e é fortemente relacionada coma etapa anterior. Aqui podemos procurar suprir as deficiências encontradas no terceiroestágio, mas deve-se também colocar em pauta outros atributos ainda não considerados.Há a necessidade de conformidade com os custos previstos, manutenção de compromissos,áreas que não devem ser preservadas tendo em vista sua importância extra-ecológica e asdiferentes preferências dentre os candidatos à preservação.
A quinta etapa é a que sai das fases de planejamento e parte para a implementaçãode unidades de conservação em si. Essa etapa é complicada especialmente pelas diversasagências, pessoas e protocolos que devem ser respeitados. E, às vezes, devido a algumimprevisto, nessa fase se fazem necessárias revisões dos planos. A implementação das áreasde conservação devem seguir dois critérios de urgência: primeiro, aqueles espaços comcaracterísticas insubstituíveis e depois áreas com grande chances de serem transformadas
1Stakeholders: termo inglês que se refere a todas as partes interessadas nos resultados de um projeto.
6
ou degradadas por ações extrativas. Aquelas áreas que possuem esses dois atributos devemter seus processos privilegiados.
A sexta etapa, e final, é a manutenção e monitoramento das reservas e provavelmentea mais longa das etapas. Dentre as tarefas desse estágio, encontramos a obrigação demanter o valor natural da área, proteger de eventuais problemas externos a ela. Duranteessa etapa deve-se periodicamente refazer todas as etapas anteriores à fim de manteratualizadas as informações sobre a reserva.
2.2 Ferramentas computacionais
Com o crescimento da necessidade de aplicar os princípios do PSC, ferramentas foramdesenvolvidas com propósito de facilitar esse o processo de seleção. Essa seção tem comoobjetivo mostrar algumas dessas ferramentas e como elas são aplicadas no contexto dePSC.
As três ferramentas apresentadas são Marxan2, C-Plan3 e TARGET [4]. As três sãoutilizadas para determinar possíveis áreas de conservação, com métodos distintos parasolucionar o problema. Marxan e C-Plan utilizam recozimento simulado (simmulatedannealing) [22, 14], mas a maneira com que esse método é simulado é diferente. Marxanprioriza custo da rede de reservas e o tamanho da fronteira de todo o sistema para alcançaros alvos de biodiversidade definidos pelo usuário. Por outro lado, o C-Plan inclui o conceitode irreaplaceabillity, definido em seguida. TARGET trabalha com um algoritmo gulosocujos detalhes são apresentados em uma seção a seguir. Nenhuma dessas ferramentasfaz uso de algoritmos multiobjetivo para a determinação de soluções, o que justifica aproposta de [31] implementada neste projeto.
2.2.1 Marxan
Marxan é a ferramenta para o planejamento sistemático de áreas de conservação maisconhecida e utilizada [34, 25, 33]. O aprimoramento do software continua em anda-mento [40]. A Figura 2.14 ilustra o funcionamento da interface gráfica Zonae Cogito como Marxan e o Geographic Information System (GIS) juntos.
Desenvolvida inicialmente por Ian Ball como parte de seu doutorado, supervisionadopor Hugh Possingham, na Universidade de Adelaide, Austrália, Marxan foi derivado deoutro programa, o SPEXAN, também desenvolvido por Ball, e segue muitos diretivas
2http://www.uq.edu.au/marxan/. Acesso 08 de julho de 2016.3http://www.edg.org.au/edg-free-resources/cplan.html. Acesso em 08 de julho de
2016.4http://www.uq.edu.au/marxan/tutorial/toc.html. Acesso 08 de julho de2016.
7
Figura 2.1: Zonae Cogito com Marxan e GIS.
adotadas por ele. Marxan tem como objetivo primário o de resolver o problema da co-bertura mínima de conjuntos. Nesse problema, temos a definição de alvos de conservaçãoe a aplicação deve selecionar unidades de planejamento que representam esses alvos como mínimo custo.
Marxan acha soluções para um problema matemático bem definido. O seu objetivo éclaro: minimizar o custo da rede de reservas e o tamanho da fronteira do sistema inteiroenquanto atinge os alvos de biodiversidade definidos pelo usuário. Uma fronteira grandedeve ser evitada por aumentar os custos com manutenção e gerenciamento, além de reduzira conectividade da rede.
A solução para que Marxan busque boas soluções pode ser descrita pela Equação2.2.1, sujeita as limitações determinadas pela equação 2.2.1, nas quais xi ∈ {0, 1} ∀ i.Na equação rij é o nível de ocorrência da característica j no local i, ci é o custo do locali, Ns é o número de locais, Nf é o número de características e Tj é o nível do alvo dacaracterística j. O valor de xi: 1 representa que o local foi escolhido para fazer parte dareserva e 0 não.
minimizeNs∑i
xici + bNs∑i
Ns∑h
xi(1 + xh)cvih
Nf∑i
xirij ≥ Tj ∀ j
O primeiro termo da fórmula acima é a penalidade associada ao custo de todos oslocais que estão presentes no sistema. O segundo é referente a configuração, ou formato,do sistema. Uma rede de reservas com um tamanho de fronteira menor, também cha-
8
mado de custo de fronteira, é um sistema mais compacto enquanto redes fragmentadase dispersas apresentam o custo de fronteira mais elevado. O valor de cvih é usualmentedeterminado pelo usuário para ser a distância física entre as fronteiras dos locais e porisso, normalmente, Marxan busca minimizá-la.
O segundo termo permite introduzir o custo, ou benefício, de adicionar um localparticular e outro local a qual o particular estaria "conectado". Esse é outro mecanismoque permite à rede de reservas ser mais que a simples soma de suas partes.
O parâmetro b é o multiplicador de fronteira que determina o custo do sistema dereserva relativo à penalidade espacial. A matriz CV é a matriz de conectividade com oselementos cvih, que reflete o custo da conexão compartilhado ao planejar unidades i e j.
Para criar soluções Marxan coloca os objetivos (Equação 2.2.1) e as limitações (Equa-ção 2.2.1) em uma única função objetivo e transforma as limitações em novos termosde penalidade. Isso significa que uma configuração do sistema que não atinge todos osobjetivos propostos ainda podem receber um valor. Essas penalidades incentivam, oudiminuem, o interesse da característica de conservação a ser considerada na configura-ção final. No final das contas, as penalidade indicam a dificuldade da atingir as metasestabelecidas para a característica de conservação específica.
Uma das aplicações mais conhecidas do Marxan foi seu uso durante o rezoneamento daGrande Barreira de Corais na Austrália em primeiro de julho de 2004. Esse foi um enormeprojeto que demorou uma década para a Reef Marine Park Authority completar, comas complexidades socio-políticas envolvidas. Marxan auxiliou no processo de expandir asáreas no parque onde a pesca era proibida, atingindo quase todas as metas de conservação.
A definição do custo das unidades de conservação foi derivada da combinação lineardos interesses de pesca comercial e recreativa, dos interesses indígenas e dos interesses daspessoas que preferiam que lugares específicos estivessem no sistema de reservas.
2.2.2 C-Plan
Matthew Watts e Robert Pressey começaram o desenvolvimento do C-Plan em 1995 [28].A ideia surgiu de um trabalhos anteriores sobre irreplaceability, ou a propriedade quealgo tem de ser insubstituível, em tradução livre. A pesquisa, na época, adicionou umadimensão nova no problema de cobertura mínima de conjuntos. O conceito de um sistemainterativo de sofware que apresenta opções espaciais para o gerenciamento de áreas deconservação então se formou.
Para elaborar o conceito de irreplacebility, discutimos agora os quatro aspectos quecompõem o problema de cobertura mínima de conjuntos:
9
• Unidades de planejamento, ou seja, os locais que serão considerados e analisadoscomo possíveis componentes da área de conservação a ser formada;
• Mapas das características de biodiversidade;
• Os objetivo de cada característica, por exemplo, número de hectares de um deter-minado tipo de vegetação conservado;
• Matriz de dados contendo a correlação entre o grau em que uma determinada ca-racterística ocorre em um local.
A irreplacebility surge como uma maneira de atribuir um grau de importância de umlocal relativo ao conjunto de locais estudado. Podemos notar que qualquer solução queé formada a partir de uma seleção de locais ainda representa apenas uma das possíveisalternativas, sem indicar a importância que os locais não escolhidos têm para atingir asmetas estabelecidas.
A irreplacebility atribuída à um local representa a probabilidade que um local tem deser necessário quando se busca atingir os objetivos estabelecidos. Esse valor oscila de 1.0a 0. Valores mais altos representam o caso que existem poucos ou nenhum local capazde substituir o local analisado enquanto valores mais próximos do limite inferior indicamque o local tem diversos outros possíveis locais que poderiam substituí-lo. Logo, locaiscujo valor de irreplacebilty é alto tem grandes chances de comporem a possível soluçãoproposta.
É interessante observar a natureza volátil que o valor da irreplacebility tem. Como éum valor relativo aos atributos definidos do problema, qualquer alteração em um dessesparâmetros tem impacto no seu valor. Se, por exemplo, o número de locais consideradosmudar, uma nova meta for estabelecida ou um objetivo for retirado, os reflexos nos valoresde cada um dos valores de irreplacebility também mudaria.
A necessidade de atualizar os valores de irreplacebility de acordo com as decisões feitassobre locais individuais despertou o interesse por um sistema de software interativo quepermitisse ao usuário explorar diferentes opções para atingir seus objetivos.
O protótipo do sistema de interação em tempo real foi desenvolvido em 1994. Em1995, o protótipo tornou-se um sistema totalmente operacional e finalmente em 1996ele se tornou o C-Plan, após vários ciclos de refinamento e expansão de funcionalidades.O direitos do C-Plan pertencem à New South Wales Department of Environment andClimate Change mas é disponibilizado como software livre.
C-Plan é amplamente utilizado na definição de áreas de conservação, especialmentena Austrália, sendo parte substancial dos esforços de preservação da região de Nova Galesdo Sul. No Brasil, o C-Plan foi utilizado como parte de um acordo do Banco Mundialcom o Governo do Estado do Goiás para dobrar a extensão das áreas de conservação. O
10
WWF-Brasil5 e a Conservação Internacional6 trabalharam para determinar novas áreasde conservação e definir as prioridades para implementação e restauração da vegetaçãonativa.
2.2.3 TARGET
TARGET [4] é um modulo do conjunto de software DIVERSITY e do BioRap toolbox.TARGET trabalha assumindo que todos os locais estudados apresentam pelo menos
um atributo de biodiversidade diferente. Esses atributos podem ser a presença de umadeterminada espécie, ou um habitat específico, entre outros mais. Cada atributo é relaci-onado a um ou mais locais.
Os atributos únicos têm a mesma importância durante o processo de seleção indepen-dente do tipo de dado ambiental que representa. TARGET trabalha com o conceito decomplementaridade, ou seja, a contribuição que um novo local reservado tem em relaçãoaos locais que já estão protegidos. Como a noção de complementaridade só é possívelde ocorrer quando já existe um conjunto de locais já reservados, TARGET trabalha ite-rativamente recalculando o fator da complementaridade cada vez que um novo local éadicionado.
Para a definição de prioridades, o TARGET dispõe de um método configurável. Oprimeiro passo desse método é a comparação do ganho complementar de um local como valor desse mesmo local para outros interesses, como agricultura ou extração. Isso é ochamado custo de oportunidade. Quando se compara o custo de oportunidade com o decomplementaridade temos o peso do trade-off. Esse peso pode ser ajustado pelo usuárioa fim de chegar à uma curva de trade-offs.
2.3 Otimização Multiobjetivo para PSC
O trabalho de Santos et al. [31] pode ser visto como um aprimoramento do trabalho deDiniz-Filho et al. [14]. O trabalho de Diniz-Filho propôs a utilização de alelos, deter-minados por análise molecular, como unidade básica de análise no processo de PSC. Oobjetivo do projeto era demostrar que é possível utilizar o PSC para determinar ótimasestratégias de conservação in situ7 e ex situ8 de uma única espécie.
5http://www.wwf.org.br/. Acesso 08 de julho de 2016.6http://www.conservation.org/global/brasil/Pages/default.aspx. Acesso 08 de ju-
lho de 2016.7Refere-se a estratégia de conservação na qual o sujeito a ser conservado, uma planta ou animal por
exemplo, é mantido no seu habitat natural.8Estratégia complementar a in situ, na qual o sujeito a ser conservado é retirado de seu local de origem
e estabelecido em uma unidade de conservação.
11
Figura 2.2: Distribuição das populações estudadas [31].
O baru (Dipteryx alata), uma planta típica do cerrado brasileiro, foi a espécie estu-dada. Foram coletados dados de 55 alelos de oito loci microsatélites9 codificados de 642árvores individuais em 25 populações dentro da área de ocupação. A Figura 2.2 mostraa localização geográfica das populações de baru coletadas e utilizadas no método.
O resultado encontrado após analisar os dados foi que eram necessárias conservar setedas populações estudadas para preservar toda a variabilidade genética do baru. A partirdessa análise, foi possível determinar quais populações eram mais ou menos importantespara a conservação do baru.
Na primeira etapa Santos et al. [31] utilizaram a mesma massa de dados e outraabordagem, os resultados que obtidos por Diniz-Filho et al. [14] A abordagem escolhidafoi a utilização de algoritmos de otimização multiobjetivo para determinar quais seriamas populações mais importantes.
O algoritmo escolhido foi o NSGA-II [13], que será detalhado no próximo capítulo.O NSGA-II é um algoritmo de otimização multiobjetivo genético e elitista. Em outraspalavras, é um algoritmo que busca a melhor relação entre vários objetivos distintos,representando seus resultados de forma análoga a que um código genético representa,
9Microsatélites são arranjos de nucleotídeos que se repetem na sequência de DNA.
12
distinguindo as melhores e piores soluções. Além disso, usa estratégias evolutivas paraguiar e avaliar sua busca.
Foram realizadas três diferentes variações de otimização. A primeira correspondeuà feita por Diniz-Filho: determinar o menor número de populações necessárias para en-globar toda a diversidade genética do baru. As funções de otimização utilizadas podemser descritas como minimizar o número de populações e maximizar o número de alelos aserem contemplados pela seleção das populações. Os resultados obtidos com essa otimi-zação corresponderam aos de Diniz-Filho et al. [14], além de indicar outras soluções nãoencontradas pelo primeiro método.
A segunda variação compreendia quatro funções de otimização. Além das duas dimen-sões da variação anterior, buscou maximizar a frequência total de alelos e a heterozigosedas populações locais ou maximizar o equilíbrio de Hardy-Weinberg10. Essas duas métri-cas também foram determinadas a partir do conjunto de dados de Diniz-Filho et al. [14].Assim, na segunda variação, temos duas análises distintas: uma maximizando a heterozi-gose e a outra maximizando o equilíbrio de Hardy-Weinberg.
O resultado dessa variação trouxe um subconjunto das soluções da primeira variação,mas com as novas características contempladas.
A terceira variação teve cinco funções de otimização: todas as apresentadas anterior-mente, de forma simultânea, isto é: maximizar alelos na seleção de populações, minimi-zar número de populações, maximizar a frequência de alelos, maximizar a heterozigosee maximizar o equilíbrio de Hardy-Weinberg. O resultado dessa variação também foium subconjunto das soluções encontradas na primeira variação, aquela com apenas doisobjetivos.
A conclusão do artigo coloca a abordagem multiobjetiva como uma alternativa à oti-mização mono-objetivo visto que a tarefa de determinar diferentes pesos para os váriosobjetivos de forma simultânea é complexa e talvez seja melhor delegar a escolha final aum especialista.
10Conceito que preconiza que, na ausência de pressão evolucionária, a frequência de alelos e genótiposdeverá se manter constante pelas gerações.
13
Capítulo 3
Otimização
Esse capítulo apresenta fundamentos de otimização. A Seção 3.1 descreve conceitos deotimização multiobjetivo. A Seção 3.2 apresenta algoritmos evolutivos de maneira geral.A Seção 3.3 detalha algoritmos genéticos, utilizados nesse trabalho. A Seção 3.4 descreveo algoritmo NSGA-II [13] que foi usado no presente trabalho e seu antecessor.
3.1 Conceitos Básicos
3.1.1 Otimização Multiobjetivo
Em problemas mais complexos ocorrem situações onde mais de um objetivo deve seralcançado. A busca pelo melhor custo-benefício é o clássico exemplo, o mais barato ou omais confortável? Ir pelo caminho mais curto ou pelo que passa por paisagens mais belas?
Esses são problemas que podem ser tratados com a aplicação de algoritmos que traba-lhem com otimização multiobjetivo. Nela os critérios de avaliação são múltiplos e muitasvezes antagônicos.
Para ilustrar o conceito de otimização multiobjetivo com critérios antagônicos vamosseguir um exemplo: escolha de uma viagem. Supondo que durante a busca por uma viagemperfeita, temos somente três objetivos a serem respeitados: minimizar o custo total daviagem, maximizar o conforto do trajeto e minimizar o tempo gasto no trajeto. Trêssoluções ótimas são encontradas, uma para cada critério. Um passagem em promoção,outra passagem em primeira classe e uma passagem cujo trajeto não faz paradas nemdesvios. Outras opções também apareceram, intermediárias a essas duas.
O questionamento que surge é: qual é a melhor alternativa? Fica claro que nenhumadas opções é indubitavelmente a melhor, cada uma das candidatas resolve o problema atécerto nível. Os três casos ótimos no entanto mostram bem como a priorização simples deapenas um critério faz, ou seja, negligencia o outro pelo menos parcialmente. As solu-
14
ções intermediárias também são questionáveis pois não atendem um objetivo de maneiracompleta.
Otimização Multiobjetivo e Otimização Mono-Objetivo
Existem discussões da natureza da otimização multiobjetivo [12]. Seriam as aplicaçõesmultiobjetivo na realidade de objetivo único? Essa pergunta vem da noção que, mesmotendo diversos objetivos, o que um algoritmo de otimização multiobjetivo faz na verdade écriar maneiras de classificar os critérios de avaliação entre as mais e as menos importantese extrair um só objetivo. Dessa maneira é possível transformar um problema multiobjetivoem um que possa ser tratado como objetivo único.
Por outro lado, não seriam os problemas de otimização mono-objetivo apenas um pro-blema multiobjetivo onde temos apenas um critério de avaliação? Nesse caso as questõesde trade-offs, ou seja, as decisões de não incluir parte de um critério para privilegiar outro,não existem pois temos apenas um objetivo a ser atingido.
O que distingue as duas abordagens é o fato de que, para o de otimização de únicoobjetivo, temos quase sempre uma resposta que é melhor que as outras, enquanto nosmultiobjetivos, todas as respostas que atendem aos requisitos listados são aceitáveis edevem ser consideradas de mesma importância.
Uma consequência das várias soluções ótimas em um algoritmo de otimização multi-objetivo é que, para aplicações mais práticas, devemos ter uma solução que se destaqueperante as demais, contrariando a premissa que todas as soluções são de igual valor.Temos duas abordagens para tratar esse problema.
A primeira das duas estratégias é determinar múltiplas soluções ótimas com diversostrade-offs para cada um dos critérios de avaliação. Em seguida, um usuário escolhe qualdas soluções encontradas atende da melhor maneira suas necessidades. A decisão deum usuário é essencial, pois existem critérios de alto nível que são difíceis de mensurare representar com precisão. Essas informações são baseadas em aspectos não técnicose dependentes da experiência que o usuário tem. A abordagem descrita é a chamadaprocedimento de otimização multiobjetivo ideal [12].
A segunda abordagem funciona invertendo a ordem da estratégia anterior. Primeiro,um vetor de preferências de alto nível é criado, atribuindo pesos distintos à objetivosdiferentes. Esse vetor será mais um critério na segunda parte, que é quando o algoritmobusca as múltiplas soluções, mas dentre elas uma é destacada devido aos pesos atribuídos.Espera-se que mudanças no vetor resultem em uma nova resposta.
15
Descrevendo um problema de otimização multiobjetivo
Em um problema de otimização multiobjetivo, temos um conjunto de funções objetivoque devem ser maximizadas ou minimizadas. Podemos definir um problema de otimizaçãomultiobjetivo genérico com a Equação 3.1, segundo o capítulo 2 [12]:
Maximize/minimize fm(x), m = 1, 2, 3, ...,Mi
sujeito a gj(x) ≥ 0, j = 1, 2, 3, ..., Ji
hk(x) = 0, k = 1, 2, 3, ..., Ki
x(L)i ≤ xi ≤ x
(U)i , i = 1, 2, 3, ..., n
(3.1)
Na Equação 3.1, x é um vetor de n variáveis de decisão: x = (x1, x2, ..., xn)T . A funçãoobjetivo fm(x) é a aquela para a qual devemos encontrar soluções. O último conjunto delimitações representam as fronteira que os valores de xi podem ter, variando do mínimox
(L)i até o máximo x(U)
i . Esses limites são conhecidos como espaço das variáveis de decisãoD.
Temos também as funções limitantes gj(x) e hk(x), que podem tanto serem tratadascomo “maior-ou-igual” ou “menor-ou-igual”, dependendo da formulação do problema.Uma solução que não satisfaz todas as (J+K) limitações e está dentro do espaço permitidopara as variáveis de decisão é uma solução impraticável. As soluções viáveis são aquelasque respeitam todas as restrições definidas, e juntas formam uma região factível (feasible),dentro do espaço permitido.
Em um problema dessa natureza, temosM funções, cada uma com requisitos próprios.Cada uma delas pode ser maximizada ou minimizada. No contexto da otimização, umafunção de maximização pode ser transformada em uma de minimização, multiplicando afunção objetivo por -1. Essa característica torna mais simples lidar com tipos misturadosde funções de otimização.
3.1.2 Frente de Pareto
A frente de Pareto [44] é um conceito usado no campo da otimização multiobjetivo, quedefine o grupo de soluções onde nenhuma solução é melhor do que as demais em todos osaspectos, quando todos os objetivos são considerados simultaneamente. Essas soluções sãoditas não-dominadas. As soluções que compõem a frente de Pareto têm a característicade indicarem soluções que são ótimas para um sub-conjunto de objetivos mas não paratodos simultaneamente.
A Figura 3.1(a) mostra a região factível para uma otimização em duas dimensões.Nesse exemplo é considerado que ambas as funções objetivo dessa otimização buscam
16
minimizar os seus valores. Os pontos A, E, I, O e U são soluções que estão na fronteiradessa região. O ponto Z e os outros pontos não nomeados são outras solução que estãono meio da região. A Figura 3.1(b) destaca as soluções que são indiferentes à Z, ou seja,não dominam nem são dominadas por ela. A Figura 3.1(c) destaca as soluções dominadaspor Z e na Figura 3.1(d) destaca as soluções que dominam Z. A Figura 3.1(e) destaca afrente de Pareto.
(a) Representação gráfica de uma re-gião factível, para otimização em duasdimensões. Cada ponto representauma solução.
(b) Áreas destacadas mostram as so-luções que são indiferentes à soluçãoZ para a minimização em duas dimen-sões.
(c) A área destacada mostra as solu-ções que a solução Z domina para aminimização em duas dimensões.
(d) A área destacada mostra as solu-ções que dominam a solução Z para aminimização em duas dimensões.
(e) A área destacada mostra a frentede Pareto para a minimização em duasdimensões.
Figura 3.1: Representações gráficas para otimização em duas dimensões, na qual busca-seminimizar ambas as dimensões.
17
3.2 Algoritmos Evolutivos
Algoritmos evolutivos podem ser vistos como uma das divisões da Computação Inspi-rada na Natureza, que propõe a observação da própria natureza para propor métodoscomputacionais [10]. A Figura 3.2 mostra como os algoritmos evolutivos se encaixamnessa hierarquia. Nessa figura sistemas imunológicos artificiais são utilizados para resol-ver otimização multiobjetivo inspirados nos sistemas imunológicos de organismos. Estaseção é baseada em Santos [30].
Figura 3.2: Como os algoritmos evolutivos pode ser classificados.
3.2.1 Definições
Algoritmos evolutivos são baseados nos princípios da Teoria da Evolução das Espécies e daGenética [2]. Charles Darwin descreve no seu livro “A origem das espécies” de 1859 comoo processo evolutivo funciona: indivíduos mais adaptados ao meio terão maior chance desobrevivência e portanto maior chance de passar adiante, para as próximas gerações, ascaracterísticas adquiridas de adaptação ao meio. Esse fenômeno, conhecido como seleçãonatural, ocorre quando uma população compete por uma quantidade limitada de recursose deve se proteger de predadores.
Outro aspecto crucial do processo evolutivo é que para se adaptar a mudanças no am-biente, os indivíduos sofrem mutações aleatórias, que podem ocorrer durante o processoreprodutivo e resultam em modificações passadas a geração seguinte. Essas mudançaspodem ser benéficas, permitindo que um indivíduo torne-se ainda mais adaptado, trans-mitindo a nova característica para seus descendentes, ou ruins, que diminuem sua chance
18
de sobrevivência. A tendência é que novas características benéficas sejam mantidas napopulação, enquanto que as características ruins sejam rapidamente eliminadas.
Os três aspectos de um algoritmo evolutivo são reprodução, mutação e seleção natural.A reprodução transmite as características de uma geração para a próxima, a mutaçãointroduz variedade nos novos indivíduos e a seleção natural elege os indivíduos maisadaptados para reproduzir.
A Figura 3.3, baseada na publicação Analyzing Evolutionary Algorithms [20] de Tho-mas Jansen, permite uma visão clara da progressão de um algoritmo evolutivo genérico.O primeiro passo, “Inicialização”, determina o tamanho da população, a quantidade denovos indivíduos gerados a cada iteração, a taxa de mutação (um valor que define a chanceque um novo indivíduo tem de sofrer mutação) e variabilidade genética. A “Seleção paraReprodução” escolhe os indivíduos mais adaptados para o espaço de acasalamento (ma-ting pool), depois aplica operandos de variabilidade em “Variação”, seleciona, se existirem,aqueles indivíduos que devem ser substituídos por novos em “Seleciona para Substitui-ção” e, caso a condição de parada seja atingida, chega ao fim. Caso contrário, uma novaiteração acontece.
3.2.2 Algoritmos Evolutivos e Otimização Multiobjetivo
Algumas características dos algoritmos evolutivos indicam que eles podem ser apropria-dos para realizar otimização multiobjetivo, e as diferentes soluções encontradas permitemconstruir a frente de Pareto. Além disso, os algoritmos evolutivos são flexíveis o sufici-ente para serem aplicados em diferentes situações. Problemas reais estão cada vez maiscomplexos, o que dificulta o uso de métodos exatos nas suas soluções.
Dentre os algoritmos evolutivos aplicados para resolver problemas de otimização mul-tiobjetivo temos:
• Vector Evaluated Genetic Algorithm [23] (VEGA): divide a população em n subpo-pulações, onde n é o número de objetivos a serem otimizados. Cada subpopulaçãoé avaliada em relação um objetivo, para depois serem combinadas em uma novapopulação utilizando operadores genéticos;
• Multi-Objective Genetic Algorithm [18] (MOGA): o critério de aptidão (fitness),definido em seguida, de um indivíduo é a sua não-dominação, ou seja, os indivíduosque não são dominados têm o maior valor de aptidão da população;
• Pareto Archived Evolution Strategy [23] (PAES): um arquivo é iterativamente cri-ado e atualizado com os indivíduos não dominados,o que permite comparar novosindivíduos com o critério de não-dominância. É uma estratégia (1+1), ou seja, umindivíduo pai gera apenas um filho por geração;
19
Figura 3.3: Fluxograma ilustrando algoritmo evolutivo simples, adaptada do Capítulo 2de Analyzing Evolutionary Algorithms [20]. As setas contínuas representam o fluxo decontrole e as setas tracejadas o fluxo de dados.
20
• Nondominated Sorting Genetic Algorithm (NSGA [35] e NSGA-II [13]): indivíduosnão-dominados em relação a toda a população são os primeiros a serem classificados,e depois retirados da população para que outros indivíduos dessa nova subpopula-ção sejam classificados como não-dominados, até que todos os indivíduos sejamclassificados. O NSGA-II é uma aprimoramento que introduz elitismo para pre-servar melhores soluções, além do operador de densidade para escolher melhoresindivíduos quando tem o mesmo nivel de não-dominância. Ambos os algoritmos sãoaprofundados na Seção 3.4;
• Strenght Pareto Evolutionary Algorithm (SPEA [45] e SPEA2 [43]): assim comono PAES, temos um arquivo com os indivíduos não-dominados encontrados até omomento. Cada indivíduo tem um valor de força, atribuído de acordo com o númerode outros indivíduos que o dominam, e um método é utilizado para obter uma melhordistribuição das soluções dominadas. O SPEA2 melhora alguns aspectos do SPEA,por utilizar também o número de indivíduos dominados por cada indivíduo, estimara densidade populacional com relação aos vizinhos e usar um método específicopara assegurar que as soluções pertencentes aos extremos da frente de Pareto nãose percam.
3.3 Algoritmos Genéticos
Algoritmos genéticos buscam simular o funcionamento do processo evolutivo, abstraindoos conceitos genéticos para utilizá-los em diversas aplicações. A ideia é modelar cadafator observando na natureza e assim tentar resolver problemas complexos.
Para algoritmos genéticos temos três características fundamentais [42]:
• Baseado em populações: um grupo de soluções é mantido como uma população,sendo cada solução um indivíduo que a compõe;
• Fitness-oriented: esse termo significa orientado à aptidão ou adaptação, em tra-dução livre. Cada indivíduo tem uma representação genética (código) e um valorde desempenho . O algoritmo pode privilegiar indivíduos com melhor desempenho,sendo essa a base para as características de conversão e otimização desse algoritmo;
• Motivado por variação: indivíduos passarão por diversas variações a fim de si-mular as mudanças sofridas por genes. Esse fator é essencial para buscas maisabrangentes no espaço de solução.
Depois de calcular os o valor de aptidão de cada solução, aqueles com valores maiselevados teriam maior chance de reproduzir e legar seus genes.
21
Figura 3.4: Gráfico da Equação 3.2.
Para melhor ilustrar a elaboração de um algoritmo evolutivo vamos seguir um exemplo.Na Figura 3.4 temos a representação da função definida pela Equação 3.2. O problema-exemplo para qual procuramos uma solução é: qual é o valor x que retorna o maior f(x)para 0 ≤ x ≤ 1 ?
f(x) = xsen(10πx)− xsen(5πx) (3.2)
O nosso problema pode ser escrito em como função com as restrições já explícitas,como mostrado pela Equação 3.3.
max f(x) =
f(x) = xsen(10πx)− xsen(5πx)
para 0 ≤ x ≤ 1(3.3)
Como podemos ver no gráfico da função, temos vários máximos locais. Dependendoda abordagem adotada existe a possibilidade da solução convergir para um máximo local.
Para aplicar o método dos algoritmos genéticos, o primeiro passo é criar a populaçãode possíveis soluções. E em seguida devemos nos ater a é definição de seus indivíduos.Para representar um indivíduo, podemos converter seu valor como um número binário.A precisão da notação depende da necessidade por respostas mais exatas. No caso daEquação 3.2 temos que 0 ≤ x ≤ 1 , ou seja, temos o intervalo entre 0 e 1 para determinarnossos indivíduos.
Se utilizarmos um bit para representar as respostas os valores representados serão o0 e o 1, se representarmos com dois bits, 0, 0.25, 0.5 e 1. Resumindo: a cada novo bitadicionado dobramos as possíveis respostas e assim aumentamos a precisão da resposta.O preço pago por maior precisão é o crescimento de gastos em memória e processamento,por esse motivo devemos mensurar bem nossos indivíduos. Para o nosso exemplo vamos
22
utilizar uma representação com 12 bits, ou seja, teremos 212 subdivisões entre 0 e 1. Aescolha de 12 bits para esse exemplo foi arbitrária mas essa escolha deve ser feita de acordocom a necessidade de cada aplicação.
O código binário é comumente usado para representação de indivíduos em algoritmosevolutivos, pois são uma boa analogia código genético a um cromossomo real. A Figura3.5 mostra a representação de um possível indivíduo, ou um cromossomo binário. Ovalor (0 ou 1) é chamado de alelo e sua posição relativa dentro do indivíduo (no topo) échamada localização ou locus. Essa representação com zeros e uns é chamada de genótipo,nome utilizado no campo da Genética, e indica os genes que compõem uma sequência doDNA, especificando as características que emergem desses genes. O fenótipo são essascaracterísticas e no caso do indivíduo da Figura 3.5 é o valor 0.85278320312. O cálculodo valor do fenótipo foi feito nesse caso considerando o genótipo (111111111111)2 comovalor máximo (1) e (000000000000)2 como o mínimo (0), cada incremento de uma unidadesoma-se 1
212 .
Figura 3.5: Representação de um indivíduo na população.
Depois de definir o tamanho da população e a maneira de representar nossos indivíduos,deve-se inicializá-los. Uma maneira simples é aleatoriamente determinar cada um dosvalores dos alelos em cada uma das suas posições para cada um dos indivíduos.
Como nosso objetivo é maximizar o valor da Equação 3.2, o valor de aptidão deverefletir essa diretiva. Portanto quanto maior o valor resultante que um indivíduo traz maisadaptado ele estará para o nosso caso. A função de aptidão é a própria Equação 3.2. Parao cromossomo da Figura 3.5 o resultado é 0.85278320312× sen(10π × 0.85278320312)−0.85278320312× sen(5π × 0.85278320312), ou seja, 0.220738600.
Após inicializarmos a população e calcularmos o valor de aptidão de cada um dosindivíduos pertencentes a ela passamos para o processo de seleção. Nesse exemplo, aquelesindivíduos cujos valores de aptidão são os mais altos devem ter maiores chances de seremselecionados. O espaço de memória onde os indivíduos selecionados vão é denominadoespaço de acasalamento. O pensamento por trás do espaço de acasalamento, assim comooutras características desse algoritmo, é o de representar a seleção natural, indivíduos maisadaptados tem maiores chances de reproduzir. Depois de definir o espaço de acasalamento,deve-se fazer a seleção entre os candidatos à reprodução.
23
Figura 3.6: Figura esquemática de um torneio binário.
Existem vários métodos diferentes de fazer a seleção dos indivíduos de uma população,a seleção por roleta, por exemplo, atribui chances maiores de seleção para indivíduos quetenham melhores valores de aptidão relativos aos demais candidatos.
Outro método que pode ser aplicado é o torneio binário, ilustrado na Figura 3.6.Temos os candidatos, no espaço de acasalamento, que são selecionados aleatoriamentepara participarem de um torneio. No caso da Figura 3.6 temos A, D, P e J comoparticipantes. Os participantes do torneio devem comparar seus valores de aptidão eaqueles com os melhores seguem para a próxima etapa. Os semifinalistas devem competirnovamente e definir o vencedor. O vencedor é o selecionado para reprodução. Na Figura3.6, o campeão foi o candidato P.
O torneio binário pode sofrer alterações para se adequar a diferentes aplicações. Onúmero de etapas de competição pode crescer ou diminuir, o critério de desempate podeser definido de diversas maneiras, um candidato pode participar de um mesmo torneio emdiferentes entradas, pode ou não haver a retirada dos participantes perdedores do espaçode acasalamento.
Alguns fatores devem ser respeitados durante a fase de seleção, como a aleatoriedadena seleção dos candidatos e um critério claro para selecionar os melhores candidatos. Amaneira como esses fatores são reforçados é versátil e deve se adaptar as necessidades doproblema. Por exemplo podemos ter um maior o menor número de etapas de avaliação,comparar os novos candidatos com um arquivo com os melhores indivíduos das etapas
24
passadas.A próxima etapa do algoritmo é a aplicação de operadores de variação a fim de gerar
novas soluções a partir das selecionadas no espaço de acasalamento. Vamos explicar duastécnicas: crossover, chamada de recombinação, e mutação [42].
A recombinação é um fenômeno que ocorre na natureza, no qual cromossomos trocamgenes de partes de DNA. O algoritmo emprega uma técnica semelhante que adota omesmo nome. Seleciona-se dentro do espaço de acasalamento dois indivíduos para seremos cromossomos pais. Essa seleção pode ser feita pareando dois elementos aleatoriamenteou determinando um critério a priori. A recombinação de ponto único troca os genesdos pais a partir de um determinado locus, muitas vezes o locus escolhido é a metade dotamanho do cromossomo [42]. Outros tipos de recombinação podem ser implementadoscomo o de dois pontos onde o espaço compreendido entre dois pontos é trocado entre doisindivíduos. Como podemos ver na Figura 3.7, a partir da posição 5 há a troca de genes.
Figura 3.7: Recombinação de ponto único.
A outra técnica de variação, a mutação, também foi criada de baseada em conceitosbiológicas. No caso natural, por algum acidente, se um gene foi mal copiado duranteuma divisão celular, por exemplo, mudanças são introduzidas em novo código genético.Para representar essa mutação pode-se fazer uso da maneira como os cromossomos foramrepresentados, ou seja, pelo código binário. A chamada bit-flip mutation, ou mutação detroca de bit em tradução livre, inverte um determinado gene de 0 para 1 e vice-versa. Aposição do gene a ser trocado poder ser determinada aleatoriamente, como pode ser vistona Figura 3.8. Podemos ressaltar que essas duas técnicas podem ser usadas em conjuntoou de maneira separada.
Depois de realizarmos as operações de variação devemos ter uma nova populaçãoderivada da população original. Comparamos os valores de aptidão dos novos indivíduoscom os da original e selecionamos o tamanho da população de indivíduos com melhor valorde adaptação. Repetimos esse processo de criar novas novas gerações até que o critério de
25
Figura 3.8: Mutação bit-flip no gene 6.
parada seja atingido, por exemplo, número máximo de gerações ou a convergência paraum conjunto de soluções.
3.4 Algoritmo NSGA-II
3.4.1 NSGA
O Nondominated Sorting Genetic Algoritm (NSGA) [35], algoritmo genético de ordenaçãonão-dominada, em tradução livre, é a primeira versão do NSGA-II utilizado nesse projeto.O funcionamento do NSGA segue a mesma linha de algoritmos genéticos simples mas édiferente na fase seleção.
Antes de ocorrer a seleção todos os indivíduos são alocados em ranks de acordo com anão-dominação desse indivíduo. Uma solução x é dominante quando comparada a outrasolução y quando todos os resultados das funções de otimização são melhores para x quepara y. Se em qualquer uma das funções, o resultado de y for melhor do que x, nenhumadas duas soluções é dominante ou dominada. No caso que de y ser superior em todos osresultados, então y domina x.
Aqueles não-dominados formam o primeiro front. Um front é um subconjunto dapopulação que agrupa indivíduos com mesmo nível de não dominância. Os indivíduos doprimeiro front recebem um valor de dummy fitness, chamado a partir de agora de aptidão-fantoche, alto. O valor de adaptação-fantoche representa o front que o indivíduo pertence,um valor arbitrário mais alto indica que o indivíduo faz parte de um dos primeiros fronts.Todos os indivíduos de uma mesmo front recebem o mesmo valor de adaptação-fantoche.
Depois de terem seu valores de adaptação-fantoche definidos há o compartilhamentoque é algo semelhante a normalização entre vizinhos. O compartilhamento busca mantera diversidade em uma população, dividindo o valor de fitness original de um indivíduoentre seus vizinhos por uma quantidade proporcional a eles.
26
Os indivíduos do front atual são ignorados após o compartilhamento e outros indiví-duos não-dominados são buscados no subconjunto restante. Repete-se o processo até nãohaverem indivíduos sem fronts.
A seleção ocorre de acordo com o valor de adaptação-fantoche. Os indivíduos comvalores mais altos tem maiores chances de serem selecionados e reproduzir. Em seguidaas operações para aumentar a diversidade são realizadas. E se for a geração máximatermina a execução, caso contrário uma nova iteração deve ocorrer. A Figura 3.9 mostrao funcionamento do NSGA.
3.4.2 Descrição do NSGA-II
O NSGA-II [13] é o aprimoramento do algoritmo genético de ordenação não-dominada ebusca solucionar alguns problemas que o seu anterior apresentava.
Uma das novas características do NSGA-II é o emprego do elitismo. Elitismo emum algoritmo evolutivo discrimina soluções em melhores e piores, dependendo de suaperformance. As soluções que apresentam resultados melhores são priorizadas duranteas outras fases do algoritmo, em especial durante a seleção para formar o espaço deacasalamento.
Para explicar os passos importantes do NSGA-II, temos que passar por três etapas:como é feita a Fast Nondominated Sorting Approach1, como preservar a diversidade dassoluções e como o laço principal é definido.
Fast Nondominated Sorting Approach
Essa abordagem busca resolver alguns problemas que o NSGA tem. O primeiro passo écalcular duas entidades para cada solução p:
• O contador de dominação np representa o número de soluções que dominam a soluçãop;
• O conjunto de soluções Sp lista todas as soluções que p domina.
Depois de calcular esses dois valores, devemos separá-las em fronts. Fronts são asestruturas responsáveis por manter o elitismo nesse algoritmo, ou seja, um front agrupasoluções em níveis hierárquicos. As melhores soluções ficam em um front dedicado, ficandoas soluções não tão boas no segundo front, e assim por diante.
O primeiro front é determinado pela soluções que tiverem np igual a zero, ou seja, nãosão dominadas por outra solução. Em seguida para cada solução p com np = 0 visitamoscada um dos membros (q) de seu Sp e diminuímos o valor do seu np pelo valor 1. Se o novo
1Abordagem de ordenação não-dominada rápida, em tradução livre.
27
Figura 3.9: Um conjunto de soluções para otimização 5D, resultante de uma execução.Visão do número de populações selecionadas, número de alelos faltantes e frequênciaalélica.
28
np de q for zero colocamos q numa nova lista Q, que formará o segundo front. Repetimosesse processo com Q, para determinar o terceiro front. Realizamos esse processo, até quetodos os indivíduos tenham seus fronts.É possível que exista o mesmo número de frontsdo que de soluções. Nesse caso, cada front contém uma única solução
Preservação de Diversidade
A maneira com que o NSGA-II mantém a variedade das soluções é composta por duasestratégias: estimação de densidade e operador de crowded-comparison 2. Para obteruma estimativa da densidade de soluções ao redor de uma solução particular, calculamosa distância média de dois pontos ao lado da solução para cada objetivo. Essa distância,chamada de crowding distance [13], é usada para estimar o perímetro do cuboide formadousando os vizinhos mais próximos como vértices. A crowding distance mostra o quãoperto uma solução está das outras soluções mais próximas à ela. Soluções muito próximassão consideradas similares entre si, enquanto soluções distantes são diferentes.
Para calcular a crowding distance, precisamos ordenar a população em ordem ascen-dente para cada uma das funções objetivo. Para cada função, as soluções nas extremidades(a de menor valor e a de maior) são as soluções de fronteira, sendo suas distâncias de-finidas como infinitas. As outras soluções intermediárias tem seus valores de distânciacalculadas como a diferença absoluta normalizada dos valores das funções das suas duassoluções adjacentes. Esse cálculo é feito em todas funções de objetivo para cada uma dassoluções. A crownding distance total é a soma de cada um dos valores de distância paracada função objetivo. As funções objetivo devem ser normalizadas, antes de calcular acrowding distance. A Equação 3.4 mostra o cálculo da crowding distance. I[i]distância é ovalor de crowding distance para a solução na posição i. I[i+ 1] e I[i− 1] são os valores decrowding distance dos vizinhos de I[i] em para o objetivo m. Finalmente, fmax
m e fminm são
os valores de fronteira para o objetivo m. Esse cálculo é feito para todos os m objetivos.O valor da crowding-distance é calculado por:
I[i]distância = I[i]distância + I[i+ 1].m− I[i− 1].m(fmax
m − fminm ) (3.4)
Quando cada solução tem seu valor de crowding distance calculado, é possível avaliaras soluções do ponto de vista de diversidade. Crowding distance pequena significa que asolução está numa área densa de soluções parecidas, enquanto um valor alto é indicativode uma solução mais diferente.
O operador de crowded-comparison (≺n) guia o processo de seleção em vários estágiosdo algoritmo para uniformemente espalhar a frente de Pareto ótima.
2Crownded-comparison: comparação de densidade, em tradução livre.
29
Para cada indivíduo i de uma população sendo irank o rank de não-dominação e idistance
o valor da crowding distance temos as comparações definidas a seguir:
i ≺n j se (irank < jrank)
OU ((irank = jrank) E (idistance > jdistance))
Dessa maneira, respeitamos o elitismo enquanto mantemos uma variedade alta. Umasolução de rank inferior não é escolhida enquanto os fronts melhores ainda não estiveremesgotados e assim soluções diferentes são priorizadas.
Laço principal
O laço principal define a iteração do algoritmo. A Figura 3.10 auxilia na compreensãodas etapas. Inicialmente uma população Pt (a população mãe) é criada aleatoriamente,é a fase 1 na figura. A população é ordenada em relação ao seu fator de não-dominação(np). Utilizando técnicas de seleção por torneio binário, recombinação e operadores demutação criamos a população descendente Q0 (população filha) de tamanho N , a fase 2.A inicialização do algoritmo é diferente das outras iterações.
Figura 3.10: As etapas do algoritmo NSGA-II, sendo as etapas indicadas de 1 a 6.
O laço principal segue a sequência a seguir:
• A população Rt = Pt ∪ Qt de tamanho 2N é formada, Rt é a população completada fase 3;
• A população Rt é ordenada de acordo com o valor de não-dominação, entre a etapa3 e 4;
30
• O conjunto de soluções do primeiro front (F1t), as não-dominadas por nenhumaoutra solução, são privilegiadas mais que qualquer outro grupo. Se F1t for menorque N todos seus membros irão formar a nova população Pt+1 (nova populaçãomãe);
• Outros membros de Pt+1 são selecionados dos outros fronts pela ordem de seu rankaté completar N membros;
• Caso existam mais membros nos fronts do que vagas em Pt+1 o operador de crowding-comparison (≺n) define os indivíduos do front que completarão as vagas remanes-centes, essa é a fase 5 na figura;
• A nova população Pt+1 de tamanho N é estabelecida;
• Voltamos a etapa 1 com Pt+1;
• Na etapa 2, Pt+1 passa por seleção, crossover e mutação para criar uma nova Qt+1
de tamanho N ;
• Segue esse laço até o número de gerações predeterminado ser atingido.
O NSGA-II pode ser acompanhado no Algoritmo 1.
31
Algoritmo 1: NSGA-IIEntrada: Número de geraçõesSaída: Pfinal
1 início2 P0 = cria-população-inicial()3 Q0 = expande-população(P0)4 t = 15 Pt = P06 Qt = Q07 para t ≤ Número de gerações faça8 Rt = Pt ∪Qt
9 F = fast-non-dominated-sort(Rt)10 Pt+1 = ∅ e i = 111 para | Pt+1 | +Fi ≤ N faça12 crowding-distance-assignment(Fi)13 Pt+1 = Pt+1 ∪ Fi
14 i = i+ 115 fim16 Ordene(Fi,≺n)17 Pt+1 = Pt+1 ∪ Fi[1 : (N− | Pt+1 |)]18 Qt+1 = expande-população(Pt+1)19 Pt = Pt+120 Qt = Qt+121 t = t+ 122 fim23 Pfinal = Pt+124 fim25 retorna Pfinal
32
Capítulo 4
NSGA-II usando Julia
Esse capítulo trata da linguagem utilizada para implementação do algoritmo e questões docódigo desenvolvido. A Seção 4.1 fala sobre a otimização em Julia, e a Seção 4.2 comentadetalhes da implementação.
4.1 Otimização usando Julia
Esta seção busca descrever mais profundamente a ferramenta utilizada neste projeto, alinguagem Julia. Em seguida, justificamos a sua escolha como a linguagem da nossaimplementação para o NSGA-II.
No início desse projeto decidido que, diferente da opção de [31] de elaborar o projeto emMATLAB®, toda implementação seria desenvolvida em linguagens não-proprietárias. Amotivação dessa escolha foi facilitar a eventual disseminação e prosseguimento do projetode PSC. Julia foi a eleita dentre as várias opções pois atende às necessidades do projeto,sendo mais informações discutidas nas seções a seguir.
4.1.1 Linguagem Julia
Julia1 é uma linguagem de programação dinâmica apresentada em 2012 [6] a fim de auxiliarno desenvolvimento de projetos de computação.
O desenvolvimento de algoritmos e análise de dados busca maior conveniência nassuas implementações e o reflexo dessa demanda é o crescimento em uso de linguagensde alto-nível de abstração. Ambientes que são caracterizados pela simplicidade de uso emaior produtividade como o MATLAB®, Octave [27], R [19] e ScyPi [37] tem ganhadodestaque para aplicações nos campos de Engenharia [1, 3], Matemática Aplicada [16, 17]
1http://julialang.org
33
e outras ciências. Entretanto, o que esses ambientes trazem em benefícios de legibilidadee capacidade de escrita, compensam na falta de desempenho.
Para resolver esse problema ocorre a adoção de arquiteturas híbridas [5], ou seja, duasou mais linguagens são combinadas para complementar os pontos fracos umas das outras.Aplicações que exigem eficiência têm a parte mais intensiva de computação escrita emlinguagens que sabidamente possuem maior eficiência, como C e Fortran [5].
Enquanto a linguagem adotada pelos desenvolvedores for mais complicada do que aempregada pelos usuários, a computação numérica será prejudicada. Surge o “problemadas duas linguagens” [5], isto é, desenvolvedores devem juntar pelo menos duas linguagenspara criar aplicações numéricas: uma excepcional para elaboração de lógica em alto nível eoutra com boa performance computacional. Essa tarefa pode ser difícil de ser completadapois, ao acoplar duas linguagens, o resultado pode não ser tão eficiente quanto se espera.
Julia é uma linguagem que se candidata como uma solução ao “problema das duas lin-guagens”, impondo que todas as funcionalidades básicas sejam passíveis de serem escritasna própria linguagem. Assim, ela busca não forçar um programador a ter que recorrer aoutras linguagens.
Além dessa característica da implementação eficiente, Julia faz uso das técnicas mo-dernas de execução de linguagens dinâmicas. O fruto desses esforços é que Julia é capaz deassociar à velocidade de execução de uma linguagem compilada estaticamente, a flexibili-dade do comportamento interativo e a produtividade que Python [38] e Ruby apresentam.O bom desempenho de Julia pode ser atribuído a três pilares: informação rica de tipos,especialização de código para tipos definidos em run-time e compilação Just-in-time (JIT)usando o framework LLVM [24].
A Figura 4.1 mostra os resultados obtidos quando um benchmark foi criado para a ava-liação e comparação de implementações feitas em Julia com relação a outras linguagensde programação. Os testes implementaram os mesmos algoritmos em diversas lingua-gens. Os algoritmos são o cálculo recursivo da sequência de Fibonacci (fib), análise deinteiros (parse_int), quicksort (quicksort), a determinação de conjunto de Mandelbrot(mandel), somatório de π (pi_sum), cálculo de estatísticas de uma matriz aleatória(ran_mat_stat) e multiplicação de matriz aleatória (rand_mat_mul).
As linguagens comparadas foram Python 2.7.1, MATLAB® versão R2011a, Octave3.4, R 2.14.2 e JavaScript V8 3.6.6.11. C++ foi a escolha da implementação de baseline,compilada em GCC 4.2.1, selecionando o melhor tempo para cada uma das otimizações.Caso métodos nativos de operações com vetores, multiplicações de matrizes e ordenaçãoexistissem, eles seriam utilizados.
Os tempos estão normalizados em relação aos resultados obtidos na implementaçãoem C++. Todos os testes foram feitos em um MacBook Pro with a 2.53GHz Intel Core
34
Figura 4.1: Gráfico do benchmark comparando Julia com outras linguagens em escalalogarítmica [6]. O eixo das ordenadas indica é o tempo decorrido em comparação àimplementação em C e o das abscissas são os diversos algoritmos implementados emvárias linguagens.
2 Duo CPU e 8GB de 1066MHz DDR3 RAM.
4.2 Detalhes da implementação do NSGA-II em Julia
A etapa da implementação do algoritmo na linguagem Julia começa com a definição decomo uma solução é representada. A abordagem escolhida foi a criação de uma classedenominada Individual.
Um Individual é composto pelos atributos genotype, fenotype, S, n, rank e crow-dingDistance:
• Genotype: vetor de inteiros cujos valores são 0 ou 1, representa o código genéticode uma solução;
• Fenotype: vetor de floats que contém os valores derivados do genotype, esses sãoos valores usados para minimizar ou maximizar algum objetivo;
• S: lista de Individuals que são dominados pela instância atual;
• n: quantidade de Individuals que dominam a instância atual;
• rank: relacionado ao front que a instância pertence;
• crowdingDistance: fator de aglomeração das soluções, em que um valor alto indicaque a solução é mais diversa que outras com valor mais baixo.
35
Para instanciar a população de tamanho pop_size foi usado um laço que, a cadainteração, gera um novo Individual. Num primeiro momento, é definido aleatoriamenteo valor do campo genotype e, a partir desse valor, o fenotype.
O valor fenotype depende do número de objetivos a serem otimizados, sendo que, paracada objetivo, um cálculo é feito. O cálculo do valor de uma das variáveis que compõem ovetor fenotype é completamente dependente da aplicação, isto é, pode ser a resposta deuma equação pré-estabelecida cuja entrada é o genotype, pode conter valores consultadosem uma matriz, pode ser combinações dessas duas estratégias ou ainda outras abordagens.É importante ressaltar que o fenotype é resultado do genotype, assim dois indivíduoscom o mesmo genotype devem apresentar o mesmo fenotype. O valor calculado éadicionado ao vetor fenotype. Na implementação, a inicialização do fenotype ocorrelogo após a do genotype na função de inicialização do Individual. Tudo o que foidescrito até o momento corresponde a linha 2 do Algoritmo 1.
Após a população inicial P - que é um vetor de Individuals - ser definida, entra-se nomain loop do NSGA-II. O número de iterações desse loop é igual ao número de gerações(generationNumber). Isso corresponde à linha 7 do Algoritmo 1.
O primeiro passo dentro do main loop é expandir a população P. A função ex-pand_population recebe P e duplica seu número de Individuals. Durante a imple-mentação foi feita decisão de expandir a população logo no começo de cada iteração, noentanto, no Algoritmo 1 temos dois momentos onde ocorre a expansão: linhas 3 e 18.Nessa etapa é introduzida mais diversidade a P. Escolhe-se aleatoriamente dois Indivi-duals de P para competir em um torneio binário de apenas uma etapa, para ser um doispais de um novo Individual. O critério para definir o vencedor do torneio é o mesmousado para definir dominância: um Individual melhor em todas as funções objetivo doque o outro ganha. Caso nenhum dos dois domine o outro, um sorteio decide o vencedor.O segundo pai é decidido da mesma maneira.
A partir de um par de pais, os clonamos para criar um par de filhos. Depois verificamosse esses dois novos filhos devem sofrer variação. A primeira possibilidade é que elesdevam passar por um crossover, a probabilidade disso ocorrer é definida pelo valor deCROSSOVER_PROBABILITY e em seguida se cada filho deve sofrer single bit mutation,cuja a chance é definida por MUTATION_PROBABILITY. Os filhos são adicionados a P.
Após expandir P, é necessário definir os Fronts. Essa etapa corresponde à linha 9 doAlgoritmo 1. Front é outra classe da implementação e é composta por um vetor de Indivi-duals. A segmentação deP em Fronts é feita pela função fast_non_dominated_sort.Nessa função, compara-se todos os Individuals que compõem P entre si para deter-minar aqueles que não são dominados por outros Individuals e assim determinar osFronts. Os valores de n e S de cada Individual são atribuídos nesse momento. A
36
fast_non_dominated_sort retorna um vetor de Fronts, no qual o primeiro fronté composto por Individuals não dominados, enquanto os Individuals do segundo sãodominados apenas pelo primeiro, os do terceiro pelos dois primeiros, e assim por dianteaté cada Individual pertença a um front.
Com os Fronts adequadamente definidos, podemos começar a criar a nova população,a newP de tamanho pop_size, o que é feito na linha 17 do Algoritmo 1. A primeiraparte é transferir todos os Fronts que caibam totalmente na newP. Neste ponto pode-mos ter duas situações: a quantidade de Individuals adicionada a newP é exatamenteigual a pop_size; ou ainda falta adicionar alguns Individuals para completá-la. Oquestionamento que surge quando faltam Individuals é como escolher os melhores deum Front, se todos eles são equivalentes da perspectiva da dominância. Para solucionaresse problema, é calculado o valor da crowdingDistance de todos os Individuals doprimeiro Front que não foram adicionados a newP (linhas 11 a 15 do Algoritmo 1). Vejaque no Algoritmo 1 esse cálculo é feito antes de selecionarmos os Individuals para anova geração, mas na implementação ele é feito no momento que se nota a necessidadede escolher entre os Individuals de um mesmo Front.
O cálculo da crowdingDistance é feito na função crowding_distance_assigned,que recebe os Individuals de um Front e determina os seus valores de crowdingDis-tance. Após ter esses valores, ordena-se Front (linha 16 do Algoritmo 1) em ordemdecrescente de crowdingDistance e adiciona-se a newP a quantidade de Individualsfaltantes para completar os pop_size necessários. Um valor de crowdingDistance altosignifica que esse candidato tem uma solução mais “distante” das outras, isto é, umaresposta mais distinta das demais. A ideia por trás disso é preservar maior diversidadede soluções.
Neste ponto, newP tem tamanho pop_size, P é definido como newP e uma novaiteração é iniciada. Repete-se esse ciclo até o número de gerações ser atingido.
Ao final, uma população P composta pelos melhores Individuals dessa execução éobtida.
Durante o desenvolvimento desse projeto, a plataforma Github2 foi utilizada comorepositório. O código-fonte e outros artefatos podem ser encontrados em https://
github.com/gsutavo/NSGA-2-in-Julia.O código-fonte está no Anexo I.
2https://github.com/
37
Capítulo 5
Estudo de Caso: Baru
Esse capítulo apresenta os dados coletados e como esses dados foram preparados a reali-zação das simulações. A Seção 5.1 apresenta a planta baru (Dipteryx alata). A Seção 5.2descreve os dados usados e como as matrizes para otimização foram construídas. Por fim,a Seção 5.3 mostra as simulações feitas a partir das matrizes e do NSGA-II implementadoem Julia.
5.1 Baru
O barueiro, que pode ser visto na Figura 5.1(a), é uma árvore frutífera que ocorre prin-cipalmente na Região Centro-Oeste do Brasil [9]. Ela é típica das matas, cerrados ecerradões e é encontrada nos estados do Mato Grosso, Mato Grosso do Sul, Goiás e Mi-nas Gerais além do Distrito Federal. Também localiza-se no Maranhão, Tocantins, Pará,Rondônia, Bahia, Piauí e São Paulo, embora com menor frequência.
O fruto do barueiro, baru (vide Figura 5.1(b)), é trabalhado para obtenção da suacastanha, que é usada para fazer óleo, farelo e manteiga [9]. Sua polpa in natura é usadapara fazer farinha, e a produção de carvão e ácidos voláteis é feita com os resíduos dosoutros processos de aproveitamento.
O grande risco que o barueiro sofre atualmente é a procura por sua madeira, que éutilizada para fabricação de carvão vegetal, instalações de cercas, indústria moveleira econstrução civil [9]. Outro problema que a preservação do baru enfrenta é o desmata-mento, que ocorre em função do avanço da fronteira agropecuária sobre o Cerrado.
Uma abordagem para preservar a planta é a criação de reservas ambientais. O processode escolha dessas áreas é complexo mas pode ser facilitado com o uso de técnicas deotimização multiobjetivo. Em [31], um método que utiliza esse tipo de otimização paraindicar populações prioritárias no momento da definição das reservas foi desenvolvido.
38
(a) Barueiro. (b) Baru com castanha visível.
Figura 5.1: Árvore e fruto do Dipteryx alata.
5.2 Dados
Nesta seção são descritos os dados do baru recebidos dos biólogos, e em seguida comoforam obtidas as matrizes que constituem as entradas do NSGA-II.
Os dados analisados foram os mesmos de [31]: o polimorfismo de 9 loci de micros-satélites (totalizando 55 alelos distintos) de 642 árvores amostradas em 25 populaçõeslocais.
A Tabela 5.1 apresenta o nome de cada uma das 25 populações locais onde as amostrasforam feitas, o número da população e o número de árvores amostradas. O número totalde indivíduos foi de 642. Para cada indivíduo, foram sequenciados 9 loci, totalizando 55alelos encontrados. A Tabela 5.2 mostra o nome de cada um dos 9 loci, o número de alelosidentificados no locus e quais são esses alelos.
A Tabela 5.3 mostra de maneira esquemática como os dados foram agregados. Aprimeira coluna representa os indivíduos analisados, 1CMT é a primeira árvore amostradana população localizada em Cocalinho-MT cujo nome é CMT , 2CMT a segunda assimsucessivamente até a 642ª, a 30CAMT, ou seja, a 30ª árvore amostrada na populaçãolocalizada em Cáceres- MT. Para cada locus foram identificados os alelos.
A Tabela 5.4 agrupa os valores de heterozigose e homozigose identificados nas popu-lações analisadas. Os conceitos de heterozigose e homozigose indicam a semelhança dos
39
Tabela 5.1: Nome das populações amostradas e quantidade de indivíduos amostrados.
Número da Nome da Local Número depopulação população de coleta indivíduos amostrados
1 CMT Cocalinho-MT 322 ABMT Água Boa-MT 323 PCO Pirenópolis-GO 324 SMS Sonora-MS 315 AMS Alcinópolis-MS 326 ATO Alvorada-TO 327 SMGO São Miguel do Araguaia-GO 328 LGO Luziânia-GO 329 ISP Icém-SP 3110 MAMG Monte Alegre de Minas-MG 3211 ENGO Estrela do Norte-GO 1212 STGO Santa Terezinha-GO 1213 AMG Arinos-MG 3214 PMG Pintópolis-MG 3215 PMS Paraíso-MS (Chapadão do Sul) 1316 PCMS Paraíso/Camapuã-MS (Camapuã) 1317 CMS Camapuã-MS 1318 IGO Indiara-GO 1319 RAMT Araguaia-MT (Barra do Garça) 2720 RAGO Araguaia-GO (Aragarças) 3721 JGO Jandaia-GO 3222 NTO Natividade-TO 1223 ARTO Arraias-TO 1524 AQMS Aquidauana- MS 3125 CAMT Cáceres- MT 30
Total 642
alelos dos indivíduos estudados. Homozigose que há uma maior percentagem de alelosidênticos enquanto heterozigose de alelos distintos. O estudo combinado desses dois fato-res nos trazem informações dos indivíduos e da população que eles fazem parte como umtodo. Uma população com valores de heterozigose e homozigose pode ser considerada emequilíbrio mas por outro lado uma população com valor de heterozigose mais alto podeser mais interessante para ser conservada pois tem representantes vários alelos.
A partir desses dados foram derivadas as matrizes que foram usadas pelo algoritmoNSGA-II.
5.2.1 Matrizes
As matrizes são resultado da manipulação dos dados brutos e representam diferentes mé-tricas que podem ser otimizadas. As explicações de como essas matrizes foram concebidas
40
Tabela 5.2: Alelos identificados para os 9 loci sequenciados.
Locus Nº alelos AlelosBm164 8 [156] [158] [165] [168]
[170] [174] [176] [178]
DaE06 3 [212] [216] [220]
DaE12 5 [216] [218] [219] [220][222]
DaE20 4 [146] [154] [156] [158]
DaE34 11 [104] [106] [108] [110][112] [114] [116] [118][120] [122] [124]
DaE41 15 [118] [120] [122] [124][126] [128] [130] [132][134] [136] [138] [142][146] [148] [150]
DaE63 3 [208] [210] [214]
DaE67 2 [170] [176]
DaE46 4 [244] [247] [250] [253]
Tabela 5.3: Representação do parcial dos dados utilizados
Loci analisadosIndivíduo Bm164 DaE06 DaE12 DaE20 DaE34 DaE41 DaE63 DaE67 DaE461CMT [158] [158] [216] [216] [220] [220] [154] [154] [110] [110] [126] [126] [208] [208] [176] [170] [253] [253]2CMT [170] [158] [216] [216] [220] [220] [154] [154] [116] [114] [126] [126] [210] [210] [176] [170] [253] [253]
. - - - - - - - - - - - - - - - - - -
. - - - - - - - - - - - - - - - - - -
. - - - - - - - - - - - - - - - - - -29CAMT [176] [174] [220] [220] [220] [218] [154] [154] [114] [110] [132] [132] [208] [208] [176] [176] [250] [250]30CAMT [156] [156] [220] [220] [220] [218] [154] [154] [114] [110] [124] [124] [208] [208] [176] [176] [250] [250]
41
Tabela 5.4: Valores de homozigose e heterozigose para cada população analisada.
População Homozigose HeterozigoseCMT 0.249 0.322ABMT 0.291 0.382PCO 0.259 0.426SMS 0.422 0.482AMS 0.349 0.476ATO 0.382 0.488SMGO 0.354 0.336LGO 0.424 0.509ISP 0.322 0.455
MAMG 0.448 0.514ENGO 0.408 0.474STGO 0.356 0.387AMG 0.293 0.31PMG 0.176 0.257PMS 0.215 0.41PCMS 0.169 0.345CMS 0.466 0.48IGO 0.35 0.453
RAMT 0.308 0.462RAGO 0.473 0.462JGO 0.638 0.599NTO 0.458 0.393ARTO 0.381 0.372AQMS 0.442 0.459CAMT 0.278 0.466
42
serão apresentadas a seguir.
Matriz A
A Matriz A representa a presença ou ausência de determinado alelo em uma populaçãolocal. Para obtê-la foi necessário primeiramente montar uma matriz intermediária X: dapresença/ausência de alelos por indivíduo.
A Matriz X foi obtida diretamente dos dados brutos, como pode ser visto na Tabela5.5. Como cada um dos indivíduos das 25 populações tem a presença (denotada pelo valor1) ou ausência (0) de determinado alelo marcada, a Matriz X é de 642× 55 e o valor naposição Xij mostra se o indivíduo da linha i teve o alelo da coluna j identificado.
Tabela 5.5: Matriz intermediária parcial X que relaciona indivíduos com os alelos identi-ficados.
Alelos identificados
Indivíduos Bm164-1
78Bm
164-1
76Bm
164-1
74Bm
164-1
70Bm
164-1
68Bm
164-1
65Bm
164-1
58Bm
164-1
56DaE
06-2
20DaE
06-2
16DaE
06-2
12DaE
12-2
22DaE
12-2
20DaE
12-2
19
... DaE
46-2
50DaE
46-2
47DaE
46-2
44
1CMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 0 01CMT 0 0 0 1 0 0 1 0 0 1 0 0 1 0 . 0 0 02CMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 0 03CMT 0 0 0 1 0 0 1 0 0 1 0 0 1 0 . 0 0 04CMT 0 0 0 1 0 1 1 0 1 1 0 0 1 0 . 0 0 05CMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 0 06CMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 1 0 07CMT 0 0 0 0 0 0 1 0 0 1 0 0 1 1 . 0 0 0
.
.
.27CAMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 1 028CAMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 1 029CAMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 0 0 030CAMT 0 0 0 0 0 0 1 0 0 1 0 0 1 0 . 1 0 0
Para determinar a Matriz A, parcialmente representada na Tabela 5.6, os valoresda Matriz X dos indivíduos da mesma população foram agregados em uma única linha.Assim, se apenas um indivíduo de uma população apresentar o alelo, a população tambémo apresentará. O resultado é a Matriz A25×55 e o valor Aij mostra se a população i temo alelo j (valor maior que zero) ou não (valor zero).
43
Tabela 5.6: Representação parcial da Matriz A - relaciona populações com os alelosidentificados. A Matriz A transposta na íntegra pode ser encontrada no Anexo II.
Alelos identificados
População Bm16
4-1
78
Bm164-1
76
Bm16
4-1
74Bm
164-1
70
Bm16
4-1
68
Bm164-1
65
Bm16
4-1
58
Bm16
4-1
56
DaE
06-2
20
DaE
06-2
16
DaE
06-2
12DaE
12-2
22
DaE
12-2
20
DaE
12-2
19
... DaE
46-2
50
DaE
46-2
47
DaE
46-2
44
CMT 2 0 0 5 0 0 31 0 4 31 0 3 32 0 . 4 8 0ABMT 0 0 1 0 0 0 32 0 0 32 0 2 32 16 . 22 19 0PCO 0 0 0 0 1 0 32 0 0 32 0 0 32 0 . 17 2 0SMS 2 19 0 0 0 1 28 2 28 20 0 0 31 0 . 29 13 0AMS 6 2 0 0 12 2 31 0 12 27 1 1 31 2 . 30 23 0ATO 2 0 0 0 0 0 32 0 14 27 0 8 32 0 . 10 22 0SMGO 1 0 0 0 0 0 32 0 15 28 0 8 30 0 . 0 11 0LGO 0 0 0 0 0 0 32 0 7 32 0 0 32 0 . 20 6 0...
NTO 0 0 0 0 0 0 12 0 1 12 0 0 12 0 . 0 7 0ARTO 0 0 0 0 0 0 15 0 3 15 0 1 15 0 . 0 13 0AQMS 9 0 0 0 0 0 31 0 30 30 0 2 31 0 . 21 16 0CAMT 2 1 2 0 0 9 2 27 30 0 0 0 28 0 . 25 5 0
Matriz B
AMatriz B relaciona a frequência dos alelos dentro das populações locais, sendo a frequên-cia dos alelos normalizada para minimizar possíveis distorções, devido a diferentes quan-tidades de árvores amostradas em cada população.
Para elaborar a Matriz B, a Matriz A foi utilizada. O primeiro passo foi determinar afrequência de cada um dos alelos dentro de uma mesma população. Para isso, os valoresde cada linha da Matriz A foram somados e em seguida o valor de cada alelo é divididopor esse somatório. Temos:
bparcialij = aij∑25j=0 aij
Esse valor no entanto ainda não foi normalizado pela quantidade de amostras de cadapopulação local. Para fazê-lo foi necessário consultar a Tabela 5.1, que nos dá a quantidadede amostras em cada uma das populações. Finalmente, dado que i é relativo ao númeroda população da Tabela 5.1 e, simultaneamente, o valor da linha da Matriz A:
bij = bparcialij
quantidade_amostras[i]
A representação parcial da Matriz B pode ser vista na Tabela 5.7.
44
Tabela 5.7: Representação parcial da Matriz B - frequência dos alelos dentro das popu-lações locais normalizada para quantidade de amostras por população. Valores multi-plicados por 104 para melhor visualização. A Matriz B transposta na íntegra pode serencontrada no Anexo III.
Alelos identificados
População Bm16
4-1
78
Bm164-1
76
Bm16
4-1
74
Bm164-1
70
Bm16
4-1
68
Bm164-1
65
Bm16
4-1
58
Bm164-1
56... D
aE46
-250
DaE
46-2
47
DaE
46-2
44
CMT 1.77 0 0 4.43 0 0 27.52 0 . 3.55 7.10 0ABMT 0 0 0.84 0 0 0 27.10 0 . 18.63 16.09 0PCO 0 0 0 0 0.90 0 28.90 0 . 15.35 1.80 0SMS 1.67 15.88 0 0 0 0.83 23.39 1.67 . 24.23 10.86 0AMS 4.82 1.61 0 0 9.64 1.60 24.90 0 . 24.10 18.47 0ATO 1.59 0 0 0 0 0 25.57 0 . 7.99 17.58 0SMGO 0.82 0 0 0 0 0 26.38 0 . 0 9.06 0LGO 0 0 0 0 0 0 27.62 0 . 17.26 5.17 0...
NTO 0 0 0 0 0 0 70.92 0 . 0 41.37 0ARTO 0 0 0 0 0 0 56.81 0 . 0 49.24 0AQMS 7.46 0 0 0 0 0 25.70 0 . 17.41 13.26 0CAMT 1.98 0.99 1.98 0 0 8.92 1.98 26.78 . 24.80 4.96 0
45
Tabela 5.8: Representação parcial da primeira matriz intermediária para elaboração daMatriz C.
Alelos identificados
Indivíduos Bm164-1
56Bm
164-1
58Bm
164-1
65Bm
164-1
68Bm
164-1
70Bm
164-1
74Bm
164-1
76Bm
164-1
78DaE
06-2
12DaE
06-2
16DaE
06-2
20DaE
12-2
16DaE
12-2
18DaE
12-2
19DaE
12-2
20
... DaE
46-2
47DaE
46-2
50DaE
46-2
53
1CMT - 2 - - - - - - - 2 - - - - 2 . - - 22CMT - 1 - - 1 - - - - 2 - - - - 2 . - - 2
.
.
.29CAMT - - - - - 1 1 - - - 2 - 1 - 1 . - 2 -30CAMT 2 - - - - - - - - - 2 - 1 - 1 . - 2 -
Matriz C
A Matriz C relaciona a heterozigose dos alelos identificados com as populações locais.Para montar a Matriz C, tomamos os dados brutos representados na Tabela 5.3. Cada
indivíduo teve seus alelos contados e anotados em uma primeira Matriz intermediária. ATabela 5.8 mostra o resultado da contagem.
Em seguida, valores do grau de heterozigose foram associados aos resultados da Tabela5.8, da seguinte forma. Para alelos homozigóticos, aqueles com contagem 2, conferiu-se ovalor -1. Para alelos ausentes 0 e para heterozigóticos (valor 1) conferiu-se -2. O resultadodessa transformação pode ser visto na Tabela 5.9.
Com os valores do grau de heterozigose já determinados, foi necessário fazer duasagregações: juntar indivíduos de mesmas populações e juntar os alelos de um mesmolocus. O processo de agregação foi análogo ao realizado para construir a matriz A, sendoque nessa só foi feita a mescla dos indivíduos em populações. O resultado final é a matrizC25x9, visto que temos 25 populações e 9 loci analisados. O valor Cij mostra o grau deheterozigose que a população de índice i tem no determinado locus j.
A Matriz C é representada pela Tabela 5.10.
Matriz D
A Matriz D relaciona cada população com seu valor do Equilíbrio de Hardy-Weinberg(EHW). Esse equilíbrio indica que, na falta de pressão evolutiva em uma população, afrequência alélica e genética permanecerá constante pelas próximas gerações.
46
Tabela 5.9: Representação parcial da segunda matriz intermediária para elaboração daMatriz C.
Alelos identificados
Indivíduos Bm164-1
56
Bm164-1
58
Bm164-1
65Bm
164-1
68
Bm164-1
70
Bm164-1
74
Bm164-1
76
Bm164-1
78DaE
06-2
12
DaE
06-2
16
DaE
06-2
20
DaE
12-2
16
DaE
12-2
18
DaE
12-2
19
DaE
12-2
20
... DaE
46-2
47
DaE
46-2
50
DaE
46-2
53
1CMT 0 -1 0 0 0 0 0 0 0 -1 0 0 0 0 -1 . 0 0 -12CMT 0 -2 0 0 -2 0 0 0 0 -1 0 0 0 0 -1 . 0 0 -1
.
.
.29CAMT 0 0 0 0 0 -2 -2 0 0 0 -1 0 -2 0 -2 . 0 -1 030CAMT -1 0 0 0 0 0 0 0 0 0 -1 0 -2 0 -2 . 0 -1 0
Tabela 5.10: Matriz C - associa os valores de heterozigose dos alelos identificados compopulações locais.
Loci analisadosPopulações Bm164 DaE06 DaE12 DaE20 DaE34 DaE41 DaE63 DaE67 DaE46
CMT -52 -41 -41 -56 -70 -44 -77 -47 -56ABMT -35 -32 -86 -71 -74 -42 -50 -50 -95PCO -35 -32 -32 -80 -68 -50 -92 -32 -41SMS -94 -82 -40 -46 -70 -81 -70 -49 -70AMS -95 -56 -71 -58 -47 -68 -50 -47 -101ATO -38 -59 -56 -67 -77 -53 -70 -83 -98SMGO -35 -65 -50 -41 -113 -83 -44 -77 -61LGO -32 -53 -32 -86 -100 -71 -64 -32 -50ISP -56 -65 -62 -35 -83 -41 -58 -32 -71
MAMG -98 -70 -65 -101 -71 -53 -77 -32 -53ENGO -27 -20 -12 -15 -27 -21 -30 -27 -33STGO -21 -17 -12 -15 -42 -21 -33 -18 -30AMG -32 -87 -32 -44 -89 -32 -56 -56 -45PMG -35 -47 -32 -32 -50 -44 -55 -47 -52PMS -16 -18 -43 -28 -24 -22 -16 -31 -17PCMS -16 -17 -34 -19 -21 -16 -16 -16 -27CMS -33 -26 -39 -24 -36 -18 -18 -12 -48IGO -13 -18 -43 -34 -19 -18 -31 -13 -37
RAMT -30 -33 -51 -63 -57 -69 -51 -27 -45RAGO -40 -46 -76 -76 -110 -91 -112 -37 -133JGO -32 -32 -92 -77 -89 -74 -113 -32 -113NTO -12 -15 -12 -24 -33 -36 -30 -12 -33ARTO -15 -24 -18 -24 -51 -18 -36 -18 -54AQMS -58 -118 -37 -76 -70 -64 -70 -37 -81CAMT -69 -30 -60 -45 -54 -76 -53 -39 -45
47
Tabela 5.11: Valores de EHW de cada população analisada.
População Valor de EHWCMT 0.293
ABMT 0.313PCO 0.645SMS 0.142AMS 0.364ATO 0.277
SMGO -0.051LGO 0.200ISP 0.413
MAMG 0.147ENGO 0.162STGO 0.087AMG 0.058PMG 0.460PMS 0.907
PCMS 1.041CMS 0.030IGO 0.294
RAMT 0.500RAGO -0.023JGO -0.061NTO -0.142
ARTO -0.024AQMS 0.038CAMT 0.676
A definição da Matriz D nesse projeto divergiu da utilizada por [31]. Nele a Matriz Drelacionava o EHW de cada locus de uma população com as probabilidades do seu testechi-quadrado, sendo uma matriz 25 por 9.
A matriz utilizada nesse projeto foi obtida de maneira diferente. A sua elaboraçãoé feita usando os valores de homozigose e heterozigose de cada população (Tabela 5.4)e calculando: F, onde Hei é o valor de heterozigose da população i e Ho é o valor dehomozigose da mesma população. A Equação 5.1 mostra esse cálculo:
F = Hei −Hoi
Hoi
(5.1)
O resultado encontrado é uma matriz 25 por 1 que relaciona cada população comseu valor de equilíbrio. Valores altos F (mais próximos de 1) indicam uma populaçãopredominantemente homozigótica enquanto valores negativo indicam maior heterozigose.
A Tabela 5.11 mostra como a matriz D relaciona EHW com as 25 populações.
48
5.3 Simulações
As simulações seguiram a mesma sequência proposta por [31], começando com dois objeti-vos a serem otimizados e em seguida aumentando-se o número de dimensões. As próximasseções explicam como cada simulação foi realizada. Todas as execuções foram feitas emuma máquina com processador Intel® Xeon® CPU E3-1231 v3 @ 3.40GHz 3.40GHz com8.00 GB RAM no Windows 10 Home 64-bit.
5.3.1 Simulação 2D
A primeira das simulações é de apenas duas dimensões. O cálculo para definição dessesdois objetivos utiliza exclusivamente a Matriz A. Esses objetivos são:
• Minimizar o número de populações locais de D. alata;
• Maximizar o número de alelos representados nesse conjunto de populações.
Em outras palavras, procura-se o menor número de populações que inclua todos alelosidentificados.
A solução é um vetor X = x1, x2, x3, ..., x24, x25 , no qual xi ∈ 0, 1 tal que xi = 1 casoa população faça parte da solução candidata e 0 caso não. Esse vetor é o genotype daimplementação. A Figura 5.2 mostra o vetor solução, onde, caso na primeira posição ovalor seja igual a 1, temos a presença da população CMT nessa resposta.
Figura 5.2: Representação do vetor solução, com índices relativos às e as respectivaspopulações locais que representam. Neste exemplo, as populações escolhidas foram CMT,PCO, AMS, ATO, MAMG, STGO, AMG, CMS e ARTO.
As funções de otimização são definidas em seguida, onde p é o número de populações(25) e n o número de alelos (55):
min
( p∑i=1
xi
),∀j ∈ {1, 2, 3, .., n}
p∑i=1
aijxi ≥ 1 (5.2)
Por questão de simplicidade, o segundo objetivo foi transformado em uma minimizaçãodo número de alelos ausentes na solução. Dessa maneira, as funções objetivo para essaotimização 2D são as funções objetivo 5.3 e 5.4.
49
min(f1(~x)) = min(quantidade_de_populações(~x)) (5.3)
min(f2(~x)) = min(alelos_faltantes(~x)) =
min(f2(~x)) = min(#Alelos− alelos_presentes(~x))(5.4)
Os parâmetros utilizados durante a execução foram os mesmos definidos por [31] parafins de comparação:
• Tamanho da população = 500;
• Probabilidade de crossover = 0,9;
• Probabilidade de ocorrer mutação = 0,05;
• Número de gerações: 10.000.
O resultado obtido é uma população de 500 indivíduos, cada um composto por umvetor resposta (genotype). O gráfico da Figura 5.3 mostra os valores do fenotypedas soluções encontradas, no eixo das ordenadas temos o "Número de alelos faltantes",que buscamos minimizar, e no das abscissas o "Número de populações selecionadas", quetambém buscamos minimizar. Cada ponto representa uma solução encontrada em umaexecução do programa.
Para complementar a visualização apresentada na Figura 5.3 um segundo gráfico, aFigura 5.4, é apresenta o mesmo conjunto de soluções, deslocando os pontos sobrepostospara representar de maneira mais clara a concentração das soluções em alguns locais. Éimportante lembrar que duas ou mais soluções com mesmos valores de fenotype não sãonecessariamente iguais.
Tabela 5.12: Resultados encontrados para nenhum alelo faltante com as sete populaçõesselecionadas.
Solução Populações Otimização 2D em Julia [14] [31]S1 1-5-17-18-19-24-25 X X XS2 1-5-15-17-19-24-25 X X XS3 1-5-17-18-19-21-25 X XS4 1-5-15-17-19-21-25 X XS5 1-5-7-17-19-21-25 XS6 1-4-5-17-19-21-25 X
A medida utilizada para averiguar se os resultados obtidos estão corretos foi compará-los com os encontrados anteriormente por [31] e por [14] para o mesmo conjunto de dados.
50
Número de populações selecionadas
-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
-5
0
5
10
15
20
25
30
35
40
45
50
55
Núm
ero
de a
lelo
s fa
ltante
s
Resultados do NSGA-II
Figura 5.3: Um conjunto de soluções para uma execução do NSGA-II para duas dimensões.
Número de populações selecionadas
-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
-5
0
5
10
15
20
25
30
35
40
45
50
55
Núm
ero
de a
lelo
s fa
ltante
s
Resultados do NSGA-II
Figura 5.4: Um conjunto de soluções para uma execução do NSGA-II para duas dimensões.Soluções que ficariam representadas de maneira sobreposta foram deslocadas para melhorvisualização.
A Tabela 5.12 compara as três abordagens em relação aos resultados para o menor número
51
de populações que englobasse todos os alelos.A Tabela 5.13 compara esses resultados em relação a irreplaceability. Irreplaceability
é a probabilidade que uma população tem de ser selecionada, ou seja, uma populaçãocom valor de irreplaceability igual a 1 deverá certamente ser selecionada. Na Tabela5.13, o valor 1 indica que a população foi escolhida na solução Si, enquanto o valor0 indica que não foi. O valor irreplaceability é calculado baseando-se no conjunto desoluções Si encontradas, respeitando as restrições impostas, nesse caso a otimização 2Drealizada. A Tabela 5.13 mostra os valores quando o número de populações selecionadas éo menor (7 populações). Se for permitida maior flexibilidade para o número de populaçõesselecionadas (8 ou 9 populações por exemplo) os valores de irreplaceability mudariam.Com apenas 7, das 25 populações, é possível ter todos os 55 alelos localizados.
Tabela 5.13: Resultados de irreplaceability para nenhum alelo faltante com as sete popu-lações selecionadas .
Irreplaceabilitynº da População Nome da população S1 S2 S3 S4 S5 S6 [14] [31] Otimização 2D em Julia
1 CMT 1 1 1 1 1 1 1.00 1.00 1.002 ABMT 0 0 0 0 0 0 0.00 0.00 0.003 PCO 0 0 0 0 0 0 0.00 0.00 0.004 SMS 0 0 0 0 0 1 0.00 0.17 0.005 AMS 1 1 1 1 1 1 1.00 1.00 1.006 ATO 0 0 0 0 0 0 0.00 0.00 0.007 SMGO 0 0 0 0 1 0 0.05 0.17 0.008 LGO 0 0 0 0 0 0 0.00 0.00 0.009 ISP 0 0 0 0 0 0 0.00 0.00 0.0010 MAMG 0 0 0 0 0 0 0.00 0.00 0.0011 ENGO 0 0 0 0 0 0 0.00 0.00 0.0012 STGO 0 0 0 0 0 0 0.00 0.00 0.0013 AMG 0 0 0 0 0 0 0.00 0.00 0.0014 PMG 0 0 0 0 0 0 0.00 0.00 0.0015 PMS 0 1 0 1 0 0 0.50 0.33 0.5016 PCMS 0 0 0 0 0 0 0.05 0.00 0.0017 CMS 1 1 1 1 1 1 0.95 1.00 1.0018 IGO 1 0 1 0 0 0 0.50 0.33 0.5019 RAMT 1 1 1 1 1 1 1.00 1.00 1.0020 RAGO 0 0 0 0 0 0 0.00 0.00 0.0021 JGO 0 0 1 1 1 1 0.57 0.67 0.5022 NTO 0 0 0 0 0 0 0.00 0.00 0.0023 ARTO 0 0 0 0 0 0 0.00 0.00 0.0024 AQMS 1 1 0 0 0 0 0.43 0.33 0.5025 CAMT 1 1 1 1 1 1 1.00 1.00 1.00
Total 7 7 7 7 7 7 -
5.3.2 Simulações 4D
Seguindo o experimento feito em [31], a segunda simulação foi composta por duas variaçõesde 4 objetivos de otimização: a 4D-heterozigose e a 4D-EHW.
52
Além das duas otimizações apresentadas na Seção 5.3.1, temos outras três funçõesobjetivo. A primeira é comum às duas variações 4D e busca maximizar a frequênciaalélica.
Para o cálculo da frequência alélica, a Matriz B foi usada:
max(f3(~x)) = max(frequência_alélica(~x)) (5.5)
A Matriz C foi usada para maximizar a heterozigose das populações locais:
max(f4(~x)) = max(heterozigose(~x)) (5.6)
A primeira variação 4D, que utiliza valor de heterozigose, busca otimizar todas asquatro funções objetivo apresentadas até agora, as funções objetivo 5.3, 5.4, 5.5 e 5.6.
A segunda variação 4D busca otimizar o EHW, em vez da heterozigose. O conjuntode populações tenderá ao equilíbrio caso a soma dos valores da população da Matriz Dsejam próximos de zero. Assim a quinta função objetivo é:
min(f5(~x)) = min(EHW (~x)) (5.7)
A segunda variação 4D busca então otimizar as funções objetivo 5.3, 5.4, 5.5 e 5.7.A Tabela 5.14 apresenta os resultados obtidos para as duas otimizações 4D, para cada
uma das variações. Alguns aspectos devem ser ressaltados: as soluções com o menornúmero de populações (7) obtidas na otimização 2D foram encontradas para ambas asotimizações 4D, o que confere com os resultados encontrados em [31]. Nota-se que, em[31], as soluções S5 e S6 também não foram identificadas para 4D.
5.3.3 Simulação 5D
A última simulação de [31] foi utilizar todas as funções objetivo definidas ao mesmo tempo,ou seja, buscou-se minimizar o número de populações selecionadas, minimizar o númerode alelos faltantes, maximizar a frequência alélica, maximizar a heterozigose do conjuntode populações selecionado e encontrar o grupo de populações mais próximo do EHW.
As funções objetivo 5.3, 5.4, 5.5, 5.6 e 5.7 definem as funções objetivo aplicadas.A Tabela 5.15 agrega as soluções encontradas durante 10 execuções para otimização
5D. Todas as soluções nessa tabela tem todos os 55 alelos presentes. As Figuras 5.5 e 5.6mostram como as soluções de uma execução podem ficar distribuídas. Ambas figuras sãodo mesmo conjunto de soluções.
As quatro soluções com apenas sete populações encontradas anteriormente foram tam-bém identificadas nessa otimização.
53
Tabela 5.14: Tabela que agrega os resultados encontrados para duas otimizações 4D com7 e 8 populações. Todos as soluções apresentadas englobam todos alelos.
Soluções Heterozigose EHW(7 populações)
1-5-15-17-19-21-25 X X1-5-15-17-19-24-25 X X1-5-17-18-19-21-25 X X1-5-17-18-19-24-25 X X
Subtotal 4 4
(8 populações)1-3-5-15-17-19-21-25 X1-3-5-15-17-19-24-25 X1-4-5-15-17-19-21-25 X1-4-5-17-18-19-24-25 X1-5-10-17-18-19-21-25 X1-5-11-17-18-19-21-25 X1-5-12-15-17-19-24-25 X X1-5-12-17-18-19-21-25 X1-5-13-15-17-19-21-25 X1-5-14-15-17-19-24-25 X1-5-14-17-18-19-21-25 X1-5-14-17-18-19-24-25 X1-5-15-16-17-19-21-25 X1-5-15-16-17-19-24-25 X1-5-15-17-18-19-21-25 X X1-5-15-17-18-19-21-25 X1-5-15-17-18-19-24-25 X1-5-15-17-19-20-21-25 X X1-5-15-17-19-20-21-25 X1-5-15-17-19-20-24-25 X1-5-15-17-19-21-22-25 X1-5-15-17-19-21-23-25 X1-5-15-17-19-21-24-25 X X1-5-15-17-19-21-24-25 X1-5-15-17-19-22-24-25 X1-5-15-17-19-22-24-25 X1-5-15-17-19-23-24-25 X1-5-15-17-19-23-24-25 X1-5-16-17-18-19-21-25 X1-5-16-17-18-19-24-25 X1-5-17-18-19-20-21-25 X1-5-17-18-19-20-24-25 X1-5-17-18-19-21-22-25 X1-5-17-18-19-21-23-25 X1-5-17-18-19-23-24-25 X1-5-8-15-17-19-21-25 X1-5-8-17-18-19-24-25 X1-5-9-15-17-19-21-25 X1-5-9-17-18-19-24-25 X
Subtotal 19 24Total 23 28
54
Tabela 5.15: Resultados das soluções para otimização 5D para 7 e 8 populações.
Soluções7 populações 1-5-15-17-19-21-25
1-5-15-17-19-24-251-5-17-18-19-21-251-5-17-18-19-24-25
Subtotal 48 populações 1-2-5-15-17-19-21-25
1-2-5-15-17-19-24-251-3-5-17-18-19-24-251-4-5-15-17-19-21-251-5-10-15-17-19-21-251-5-10-17-18-19-24-251-5-11-15-17-19-21-251-5-12-15-17-19-24-251-5-13-15-17-19-21-251-5-15-16-17-19-21-251-5-15-16-17-19-24-251-5-15-17-18-19-24-251-5-15-17-19-21-22-251-5-16-17-18-19-24-251-5-17-18-19-21-22-251-5-17-18-19-23-24-251-5-6-17-18-19-24-251-5-8-17-18-19-21-251-5-9-17-18-19-21-251-5-9-17-18-19-24-25
Subtotal 20Total 24
55
Figura 5.5: Um conjunto de soluções para otimização 5D, resultante de uma execução.Visão do número de populações selecionadas, número de alelos faltantes e frequênciaalélica.
Figura 5.6: Um conjunto de soluções para otimização 5D, resultante de uma execução.Visão do número de populações selecionadas, valor de heterozigose e EHW.
56
Capítulo 6
Estudo de Caso: MamíferosTerrestres
Esse capítulo mostra a aplicação do método implementado para dados terrestres coleta-dos no mundo todo. A Seção 6.1 descreve os dados utilizados. A Seção 6.2 discute osresultados para o Distrito Federal e a América do Sul.
6.1 Dados dos mamíferos terrestres
Os dados dos mamíferos terrestres desse estudo de caso foram utilizados em [15]. Nesseartigo foram identificadas áreas de conservação ambiental a partir da presença/ausênciade mamíferos terrestres e projeções para o uso em agricultura no futuro.
Os dados de presença e ausência de mamíferos terrestres em todo o mundo estãoagregados em uma matriz M, representada parcialmente na Tabela 6.1. A matriz M temdimensões 87.516 (pares de coordenadas) por 5.216 (mamíferos terrestres). A posiçãoMij
mostra se o mamífero da posição j está presente, com valor 1, na área da posição i, ounão, com valor 0.
6.2 Presença e ausência de mamíferos terrestres
6.2.1 Distrito Federal e Entorno
A primeira análise dos dados da matriz M foi de uma submatriz, MDF , que engloba aárea do Distrito Federal (DF) e algumas regiões do entorno. Essa primeira análise foi feitacom um número pequeno de regiões para familiarização com o novo conjunto de dados.O método para a preparação de dados utilizado para esse caso inicial foi depois aplicadopara uma região maior.
57
Tabela 6.1: Representação parcial da matriz de presença/ausência de mamíferos terrestresem determinadas regiões geográficas. As coordenadas geográficas são dadas em notaçãodecimal.
Coordenadas Mamíferos terrestres
Long
itude
(x)
Latit
ude(y)
Abditomys
latid
ens
Abeomelom
yssevia
Abrawa
yaom
ysruschii
Abrocomabenn
ettii
Abrocomaboliv
iensis
Abrocomabudini
Abrocomacine
rea
... Zygodo
ntom
ysbrun
neus
Zygogeom
ystricho
pus
Zyzomys
argurus
Zyzomys
maini
Zyzomys
palatalis
Zyzomys
pedu
nculatus
Zyzomys
woodwa
rdi
-179,75 19,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0-179,75 18,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0-179,75 17,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0-179,75 16,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0-179,75 16,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
... ...
179,25 82,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0179,25 82,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0179,25 83,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0179,25 83,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
58
Para a elaboração da MDF foram necessários dois recortes na matriz M. O primeirofoi o de extrair somente as coordenadas da área escolhida que abrange do DF. Assim, aregião delimitada pelas coordenadas (-17,25, -47,25), (-14,75, -47,25), (-17,25, -48,75) e(-14,75, -48,75) foi recortada de M.
Após o primeiro recorte, obtivemos uma submatriz de 24 por 5.216, visto que a áreadelimitada pelas quatro coordenadas se traduz em 24 linhas da matriz M e temos osvalores de presença e ausência de todos os mamíferos terrestres contemplados pelos dados.A segunda parte da definição de MDF foi retirar todas as colunas da matriz intermediáriacuja soma de seus elementos seja igual a zero, ou seja, retirar os dados dos mamíferos quenão são encontrados no DF e Entorno.
O resultado final é a matriz MDF de 24 por 169 cuja posição MDF ij nos diz se na áreai temos o animal j. A Tabela 6.2 ilustra parcialmente como a matriz MDF é.
Tabela 6.2: Tabela parcial de valores de presença e ausência de mamíferos terrestres decada par de coordenadas do DF e Entorno.
Coordenadas Mamíferos terrestres
Posiç
ão
Long
itude
(x)
Latit
ude(y)
Akodo
nlin
dberghi
Aloua
ttacaraya
Anou
racaud
ifer
Anou
rageoff
royi
Artib
eusgn
omus
Artib
euslituratus
Artib
eusobscurus
... Thylamys
karimii
Thylamys
velutin
us
Tolypeutes
tricinctus
Tona
tiabidens
Tracho
pscirrho
sus
Uroderm
abilobatum
Uroderm
amagnirostrum
1 -47,25 -15,75 0 1 1 1 1 1 1 ... 0 1 0 1 1 1 02 -47,25 -15,25 0 1 1 1 1 1 1 ... 0 1 0 1 1 1 03 -47,25 -14,75 0 1 1 1 1 1 1 ... 0 1 1 1 1 1 04 -47.25 -14,25 0 1 1 1 1 1 1 ... 0 1 1 1 1 1 0
... ...
21 -48,75 -16,25 0 1 1 1 1 1 1 ... 0 1 0 1 1 1 022 -48,75 -15,75 0 1 1 1 1 1 1 ... 0 1 1 0 1 1 023 -48,75 -15,25 0 1 1 1 1 1 1 ... 1 1 1 0 1 1 024 -48,75 -14,25 0 1 1 1 1 1 1 ... 1 1 1 0 1 1 1
A Figura 6.1(a) mostra uma imagem aproximada da área analisada. Cada pontodestacado pela Figura 6.1(a) representa uma área de 0,5 graus de latitude e 0,5 de latitude,como ilustrado pela Figura 6.1(c). A Figura 6.1(b) mostra os índices de posição utilizadospara referênciar cada posição.
Para analisar os dados do DF e entorno, foi feita uma otimização em duas dimen-sões: minimizar o número de áreas em uma solução e minimizar o número de espécies demamíferos terrestres faltantes nas regiões selecionadas.
59
(a) Regiões do DF e entorno analisados.
(b) Mapeamentodas posiçõesde cada par decoordenadas.
(c) Legenda dos pontos nomapa. Cada ponto fica nocentro de uma área definidapor quatro coordenadas queformam retângulos de lados0,5 graus.
Figura 6.1: Área analisada durante o experimento no DF e Entorno.
60
Os parâmetros foram os mesmos definidos por [31] e utilizados nas otimizações dobaru. A motivação dessa escolha foi verificar a aplicabilidade do método definido por [31]em outro problema. Todas as execuções foram feitas em processador Intel® Xeon® CPUE3-1231 v3 @ 3.40GHz 3.40GHz com 8.00 GB RAM no Windows 10 Home 64-bit. Osparâmetros são:
• Tamanho da população = 500;
• Probabilidade de recombinação = 0,9;
• Probabilidade de ocorrer mutação = 0,05;
• Número de gerações: 10.000.
A Figura 6.2 mostra como o conjunto de respostas para cada execução do NSGA-IIpara essa aplicação. Cada ponto representa uma solução. A Tabela 6.3 agrega as soluçõesencontradas durante as execuções com o menor número de populações que representemtodos os mamíferos da região toda. O valor mínimo encontrado foi de 5 regiões, maso número de possibilidades para escolha cresce consideravelmente quando 6 regiões sãopermitidas.
Tabela 6.3: Conjunto das soluções encontradas para otimização 2D na área definida parao DF e entorno com 5 e 6 populações selecionadas. Em todos os resultados, todos osmamíferos são incluídos.
Soluções5 populações 1-6-10-13-24
1-6-10-14-241-6-10-19-241-6-10-20-241-6-10-21-24
Subtotal 56 populações 1-2-6-10-20-24
1-3-6-10-13-241-4-6-10-19-241-4-6-16-19-241-5-6-10-14-241-5-6-13-16-241-6-10-11-13-241-6-10-12-13-24
Continua na próxima página
61
Tabela 6.3 – continuando da última página1-6-10-12-20-241-6-10-12-21-241-6-10-13-14-241-6-10-13-15-241-6-10-13-17-241-6-10-13-18-241-6-10-13-19-241-6-10-13-21-241-6-10-13-23-241-6-10-14-17-241-6-10-14-19-241-6-10-14-20-241-6-10-14-23-241-6-10-16-21-241-6-10-17-19-241-6-10-17-20-241-6-10-18-19-241-6-10-19-20-241-6-10-19-21-241-6-10-20-21-241-6-10-20-23-241-6-7-10-13-241-6-7-10-19-241-6-7-10-20-241-6-8-10-14-241-6-9-10-13-241-6-9-10-19-241-6-9-10-20-24
Subtotal 36Total 41
6.2.2 América do Sul
A segunda etapa foi aplicar a mesma abordagem utilizada no DF e entorno em uma áreamaior, sendo a América do Sul a região selecionada. A escolha da América do Sul foi feitapelos motivos que o tamanho dos dados analisados nessa área é consideravelmente maior
62
Número de áreas selecionadas
-1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
-505
101520253035404550556065707580859095
100105110115120125130135140145150155160165170
Núm
ero
de e
spec
ies
falta
ntes
Figura 6.2: As soluções encontradas para otimização 2D no DF e entorno.
63
daqueles anteriormente estudados, tanto para o DF quanto para o baru, e, ao mesmotempo, ainda é pequeno quando comparado com os dados do mundo inteiro.
Assim como feito no caso do DF, o primeiro passo foi recortar a área de interesseda matriz M . As coordenadas extremas da América do Sul foram definidas a partir dopar longitude e latitude mínimas (-81,25, -56,25) e da longitude e latitude máximas (-34,75, 12,25)1. A Figura 6.3 mostra o recorte feito. Somente as áreas terrestres da regiãoselecionada foi recortada.
Figura 6.3: A América do Sul, com quatro pontos extremos em destaque, e área recortada.
Após escolher região, a segunda parte da preparação da matriz foi retirar os mamíferosque não estão presentes na América do Sul. O resultado final é a matrizMAS, parcialmentedescrita na Tabela 6.4.
Os parâmetros utilizados nesse experimento foram os mesmos que durante a otimiza-ção do DF e Entorno, salvo o número de gerações. Apenas 100 gerações foram executadas,devido ao elevado custo computacional. De uma população de 500 indivíduos, 376 repre-sentaram todos os mamíferos terrestres. A Figura 6.4 mostra a solução com o menornúmero de regiões selecionadas (3.172). Em [15], busca-se determinar áreas de conserva-ção que abranjam 17% da área terrestre, mas conservando todos os mamíferos terrestres
1Coordenadas definidas em https://worldmap.harvard.edu/maps/462.
64
cujos dados foram coletados. No nosso experimento, 46,45% da área terrestre total daAmérica do Sul foi selecionada representando todos mamíferos. Porém, esperamos que,para um maior número de gerações, esse percentual decresça.
Tabela 6.4: Tabela parcial da MAS com valores de presença e ausência de mamíferosterrestres de cada par de coordenadas da América do Sul.
Coordenadas Mamíferos terrestres
Posiç
ão
Long
itude
(x)
Latit
ude(y)
Abditomys
latid
ens
Abeomelom
yssevia
Abrawa
yaom
ysruschii
Abrocomabenn
ettii
Abrocomaboliv
iensis
Abrocomabudini
Abrocomacine
rea
... Zygodo
ntom
ysbrun
neus
Zygogeom
ystricho
pus
Zyzomys
argurus
Zyzomys
maini
Zyzomys
palatalis
Zyzomys
pedu
nculatus
Zyzomys
woodwa
rdi
1 -81,25 -6,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 02 -81,25 -5,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 03 -81,25 -5,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 04 -81,25 -4,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
... ...
6826 -34,75 -7,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 06827 -34,75 -7,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 06828 -34,75 -6,75 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 06829 -34,75 -6,25 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
65
Figura 6.4: Regiões a serem preservadas na América do Sul, com coordenadas de conser-vação destacadas.
66
Capítulo 7
Conclusões e Trabalhos Futuros
Nesse trabalho, foi implementado e analisado o método de otimização multiobjetivo ba-seado em NSGA-II para identificação de áreas de conservação proposto, em [31]. A im-plementação do algoritmo NSGA-II [13] foi desenvolvida em Julia [6] e foram realizadosdois estudos de caso.
O primeiro estudo de caso foi realizar o mesmo experimento proposto por [31], utili-zando os alelos da Dipteryx alata, o baru, planta típica do Cerrado Brasileiro. O objetivofoi assegurar a corretude da nova implementação desenvolvida para esse projeto. Os re-sultados encontrados foram consistentes com aqueles obtidos por [31] e anteriormente por[14].
O segundo estudo de caso foi aplicar o mesmo método para dados fornecidos por[15]. Esse novo conjunto de dados foi consideravelmente maior e reportava a presençade mamíferos terrestres em diversas regiões do globo terrestre. A primeira otimizaçãofoi encontrar o menor número de áreas que representasse todos os mamíferos no DistritoFederal e entorno. A segunda otimização foi representar todos os mamíferos no menorconjunto de áreas, mas em toda a América do Sul.
O próximo passo é analisar os dados obtidos até agora e depois expandir para dadosdo mundo todo. Além disso, deve-se incluir novas dimensões do problema para que osnovos resultados reflitam de melhor maneira os trade-offs relacionados à definição dereservas ambientais que existem no mundo real, em especial as informações de projeçãode áreas agrícolas seguindo o que foi proposto em [15]. Temos a intenção de aprimoraro próprio método de otimização, comparando uma implementação paralela do NSGA-II [41] com a desenvolvida neste trabalho. Pretendemos ainda investigar paralelismo ecomputação distribuída usando Julia. Por fim, gostaríamos de construir uma plataformade fácil acesso, que auxilie tomadores de decisão a definir áreas de conservação com baseem critérios científicos.
67
Referências
[1] A. Azemi e L.L. Pauley. Teaching the introductory computer programming course forengineers using MATLAB. In Frontiers in Education Conference, 2008. FIE 2008.38th Annual, pages T3B–1. IEEE, 2008. 33
[2] T. Back. Evolutionary algorithms in theory and practice: evolution strategies, evolu-tionary programming, genetic algorithms. Oxford university press, 1996. 18
[3] S.I. Barry e T. Webb. Multi-disciplinary approach to teaching numerical methods toengineers using MATLAB. Anziam Journal, 47:216–230, 2006. 33
[4] D. Barton, G. Rusch, J.O. Gjershaug, D.P. Faith, e L. Paniagua. TARGET as atool for prioritising biodiversity conservation payments on private land-a sensitivityanalysis. NIVA-rapport, 4856, 2004. 7, 11
[5] J. Bezanson, A. Edelman, S. Karpinski, e V.B. Shah. Julia: A fresh approach tonumerical computin. arXiv preprint arXiv:1411.1607, 2014. 34
[6] J. Bezanson, S. Karpinski, V.B. Shah, e A. Edelman. Julia: A fast dynamic languagefor technical computing. arXiv preprint arXiv:1209.5145, 2012. ix, 33, 35, 67
[7] C.A. Bianchi e S.M. Haig. Deforestation trends of tropical dry forests in centralbrazil. Biotropica, 45(3):395–400, 2013. 1
[8] P.M. Brando, M.T. Coe, R. DeFries, e A.A. Azevedo. Ecology, economy and ma-nagement of an agroindustrial frontier landscape in the southeast amazon. Phi-losophical Transactions of the Royal Society of London B: Biological Sciences,368(1619):20120152, 2013. 1
[9] L.R. Carrazza e J.C.C. D’Ávila. Manual tecnológico de aproveitamento integral dofruto do baru (dipteryx alata). Technical report, Instituto Sociedade, População eNatureza (ISPN), 2010. 38
[10] L.N. de Castro. Fundamentals of natural computing: an overview. Physics of LifeReviews, 4(1):1–36, 2007. 18
[11] G. de Oliveira, M.B. Araújo, T.F. Rangel, D. Alagador, e J.A.F Diniz-Filho. Con-serving the brazilian semiarid (caatinga) biome under climate change. Biodiversityand Conservation, 21(11):2913–2926, 2012. 6
[12] K. Deb e D. Kalyanmoy. Multi-Objective Optimization Using Evolutionary Algo-rithms. John Wiley & Sons, Inc., New York, NY, USA, 2001. 15, 16
68
[13] K. Deb, A. Pratap, S. Agarwal, e T.A.M.T. Meyarivan. A fast and elitist multiobjec-tive genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation,6(2):182–197, 2002. 12, 14, 21, 27, 29, 67
[14] J.A.F. Diniz-Filho, D.B. Melo, G. de Oliveira, R. G. Collevatti, T.N. Soares, J.C. Na-bout, J. de Souza Lima, R. Dobrovolski, L. J. Chaves, R.V. Naves, et al. Planning foroptimal conservation of geographical genetic variability within species. ConservationGenetics, 13(4):1085–1093, 2012. 7, 11, 12, 13, 50, 52, 67
[15] R. Dobrovolski, R. Loyola, G.A.B. Fonseca, J.A.F. Diniz-Filho, e M.B. Araújo. Glo-balizing conservation efforts to save species and enhance food production. BioScience,64(6):539–545, 2014. 57, 64, 67
[16] J.J. Faraway. Linear models with R. CRC Press, 2014. 33
[17] B. Fischer, J. Schumann, e T. Pressburger. Generating Data Analysis Programs fromStatistical Models. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000. 33
[18] C.M. Fonseca e P.J. Fleming. Genetic algorithms for multiobjective optimization:Formulation, discussion and generalization. Proceedings of the 5th International Con-ference on Genetic Algorithms,, 63:416–423, 1993. 19
[19] R. Ihaka e R. Gentleman. R: a language for data analysis and graphics. Journal ofcomputational and graphical statistics, 5(3):299–314, 1996. 33
[20] T. Jansen. Analyzing evolutionary algorithms: The computer science perspective.Springer Science & Business Media, 2013. ix, 19, 20
[21] A.K. Jorgenson e T. Dietz. Economic growth does not reduce the ecological intensityof human well-being. Sustainability Science, 10(1):149–156, 2015. 1
[22] S. Kirkpatrick. Optimization by simulated annealing: Quantitative studies. Journalof statistical physics, 34(5-6):975–986, 1984. 7
[23] J.D. Knowles e D.W. Corne. Approximating the nondominated front using the paretoarchived evolution strategy. Evolutionary computation, 8(2):149–172, 2000. 19
[24] C. Lattner e V. Adve. Llvm: A compilation framework for lifelong program analy-sis & transformation. In Code Generation and Optimization, 2004. InternationalSymposium on Code Generation and Optimization, pages 75–86. IEEE, 2004. 34
[25] R. Lourival, H. McCallum, G. Grigg, C. Arcangelo, R. Machado, e H. Possingham.A systematic evaluation of the conservation plans for the pantanal wetland in Brazil.Wetlands, 29(4):1189–1201, 2009. 7
[26] C.R. Margules e R.L. Pressey. Systematic conservation planning. Nature,405(6783):243–253, 2000. 5
[27] M. Murphy. Octave: A free, high-level language for mathematics. Linux journal,1997(39es):8, 1997. 33
69
[28] R.L. Pressey, M.E. Watts, T.W. Barrett, e M.J. Ridges. The C-plan conservationplanning system: origins, applications, and possible futures. Spatial conservationprioritization: quantitative methods and computational tools, pages 211–34, 2009. 9
[29] A.A. Rogers, J.A. Cleland, e M.P. Burton. The inclusion of non-market values in sys-tematic conservation planning to enhance policy relevance. Biological conservation,162:65–75, 2013. 6
[30] S.S. Santos. Otimização multiobjetivo aplicada ao planejamento sistemático de con-servação para espécies de plantas do cerrado brasileiro. Tese (Doutorado), Universi-dade de Brasília, 2015. 18
[31] S.S. Santos, M.E.M.T. Walter, A.C.P.L.F. Carvalho, T.N. Soares, M.P.C. Telles, R.D.Loyola, e J.A.F. Diniz-Filho. Multi-objective optimization in systematic conservationplanning and the representation of genetic variability among populations. Geneticsand Molecular Research, 14(2):6744–6761, 2015. ix, 2, 4, 7, 11, 12, 33, 38, 39, 48, 49,50, 52, 53, 61, 67
[32] S.S. Santos, M.E.M.T. Walter, A.C.P.L.F. de Carvalho, T.N. Soares, M.P.C. Telles,R.D. Loyola, e J.A.F. Diniz-Filho. Multi-objective optimization for plant germplasmcollection conservation of genetic resources based on molecular variability. Tree Ge-netics & Genomes, 11(2):1–10, 2015. 6
[33] J. de A. Silva, R.B. Machado, A.A. Azevedo, G.M. Drumond, R.L. Fonseca, M.F.Goulart, E.A.M. Junior, C.S. Martins, e M.B. Ramos-Nero. Identificação de áreasinsubstituíveis para conservação da Cadeia do Espinhaço, estados de Minas Gerais eBahia, Brasil. Megadiversidade, 4:272–309, 2008. 7
[34] R.J. Smith, P.D. Eastwood, Y. Ota, e S.I. Rogers. Developing best practice for usingmarxan to locate marine protected areas in european waters. ICES Journal of MarineScience: Journal du Conseil, 66(1):188–194, 2009. 7
[35] N. Srinivas e K. Deb. Muiltiobjective optimization using nondominated sorting ingenetic algorithms. Evolutionary computation, 2(3):221–248, 1994. 21, 26
[36] L. Starke, E. Assadourian, e M. Renner. State of the World 2012: Moving TowardSustainable Prosperity. Island Press/Center for Resource Economics, 2012. 1
[37] S. Van Der Walt, S.C. Colbert, e G. Varoquaux. The NumPy array: a structure forefficient numerical computation. Computing in Science & Engineering, 13(2):22–30,2011. 33
[38] G. van Rossum e J. de Boer. Interactively testing remote servers using the pythonprogramming language. CWi Quarterly, 4(4):283–303, 1991. 34
[39] S. Veloz, L. Salas, B. Altman, J. Alexander, D. Jongsomjit, N. Elliott, e G. Bal-lard. Improving effectiveness of systematic conservation planning with density data.Conservation Biology, 29(4):1217–1227, 2015. 6
70
[40] M.E. Watts, I.R. Ball, R.S. Stewart, C.J. Klein, K. Wilson, C. Steinback, R. Lourival,L. Kircher, e H.P. Possingham. Marxan with zones: software for optimal conservationbased land-and sea-use zoning. Environmental Modelling & Software, 24(12):1513–1521, 2009. 7
[41] S. Yijie e S. Gongzhang. Improved NSGA-II multi-objective genetic algorithm basedon hybridization-encouraged mechanism. Chinese Journal of Aeronautics, 21(6):540–549, 2008. 67
[42] X. Yu e M. Gen. Introduction to evolutionary algorithms. Springer Science & BusinessMedia, 2010. 21, 25
[43] E. Zitzler, M. Laumanns, L. Thiele, et al. SPEA2: Improving the strength paretoevolutionary algorithm. Eurogen, 3242(103):95–100, 2001. 21
[44] E. Zitzler e L. Thiele. An evolutionary algorithm for multiobjective optimization:The strength pareto approach. Technical Report 54, Computer Engineering andCommunication Networks Lab (TIK), Swiss Federal Institute of Technology (ETH),1998. 16
[45] E. Zitzler e L. Thiele. Multiobjective evolutionary algorithms: a comparative casestudy and the strength pareto approach. IEEE transactions on Evolutionary Com-putation, 3(4):257–271, 1999. 21
71
Anexo I
Código-fonte do NSGA-II em Julia
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Individual Type
##
"""
Individual Type represents a possible solution
"""
type Individual
# represents this solution’s genetic code: an array
# of integers 0 or 1
genotype::Array{Int8,1}
# represents the features that emerge from this
# solution’s genotype
fenotype::Array{Float32}
# array of Individuals dominated by this one,
# initialized empty
S::Array{Individual,1}
# number of Individuals that dominate this one,
# initialized 0
n::Int
# Individual’s rank
rank::Int
# Crowding distance value
72
crowdingDistance::Float32
# Individual’s constructor, receives an integer and
# initializes a genotype with random values
function Individual(size::Int)
genotype = initGene(size)
fenotype = initFenotype(genotype)
new(genotype,fenotype ,[],0,0,0)
end
# Individual’s constructor, receives an array of
# integers and initializes a genotype with it
function Individual(gene::Array{Int8})
genotype = copy(gene)
fenotype = initFenotype(genotype)
new(genotype,fenotype,[],0,0,0)
end
end
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Guilherme N. Ramos ([email protected])
# NSGA-II
##
include("Individual.jl")
include("Front.jl")
include("variation.jl")
include("initialization.jl")
include("display.jl")
include("utils.jl")
"""
73
Assign some constants
"""
data_matrixA = readcsv("../data/matrizA.csv")
data_matrixB = readcsv("../data/matrizB_normal.csv")
data_matrixC = readcsv("../data/matrizC_25by9.csv")
data_matrixD = readcsv("../data/matrizD.csv")
geneSize = 25
pop_size = 500
CROSSOVER_PROBABILITY = 0.9
MUTATION_PROBABILITY = 0.05
generationNumber = 10000
"""
Runs NSGA-II
"""
function nsga2()
newP::Array{Individual} = []
# Population of size pop_size is created
P = initPopulation(pop_size)
for k = 1:generationNumber
if k == 1 || k == 2500 || k ==5000 || k==7500
println("Generation ", k)
end
# Initial populations expands to double
# its initial size
expand_population(P)
newP = []
# Set ranks and fronts
F = fast_non_dominated_sort(P)
i = 1
74
while length(newP) + length(F[i].individuals) <= pop_size
append!(newP, F[i].individuals)
i = i + 1
end
if length(newP) < pop_size
# Set the crowding distance of a front’s
crowding_distance_assigned(F[i].individuals) individuals
# Sort the front based on crownding distance
sort!(F[i].individuals,
lt = (x,y)-> x.crowdingDistance < y.crowdingDistance,
rev = true)
individualsNeeded = pop_size - length(newP)
for j in 1:individualsNeeded
if F[i].individuals[j].crowdingDistance > 0
newP = push!(newP, F[i].individuals[j])
else
append!(newP,
random( individualsNeeded - j,F[i].individuals[j:end]))
break
end
end
end
P = copy(newP)
end
println("-----------------------------------")
println("Results of interest:")
println("-----------------------------------")
sort!(P, lt = (x,y)-> x.crowdingDistance < y.crowdingDistance,
rev = true)
# Print file with individuals that attend
# some predefined parameters
debug_printFile(P)
75
# Create a intermediate file to plot graphs
createInterFile(P)
end
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Front Type
##
"""
Front Type represents a sub-group a population.
A population may be divided as fronts, the first front
englobes individuals that dominate all others.
The second front, the individuals that are only dominated
by the ones in the first front and that dominate all
other individuals.
The third is dominated by the first and second and so on.
There can be as many fronts as individuals.
"""
type Front
individuals::Array{Individual,1}
function Front(pop::Array{Individual})
individuals = copy(pop)
new(individuals)
end
function add(new::Individual)
individuals = push!(individuals,new)
end
end
########################################################
76
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Initialization functions
##
"""
Initializes a gene, receives the gene size (size)
and returns an array with random values 0 or 1
"""
function initGene(size::Int)
array::Array{Int8,1} = []
array = rand(0:1,size)
return array
end
"""
Initializes the fenotype. receives an array of integers
(should be the genotype) and calculates the features
based on the array received
Returns an array of features
"""
function initFenotype(entry::Array)
objetiveValue = 0
# Defines how many different alelles are expected
sizeAlelleArray = size(data_matrixA)[2]
exit::Array{Float32} = []
auxArray::Array{Int32} = fill(0,sizeAlelleArray)
"""
First objetive: number of 1’s in genotype
"""
objetiveValue = sum(entry)
77
exit = push!(exit, objetiveValue)
"""
Second objetive: minimize missing alleles
"""
for i = 1:length(entry)
if entry[i] == 1
for j = 1:sizeAlelleArray
if data_matrixA[i,j] > 0
auxArray[j] = 1
end
end
end
end
objetiveValue = 0
x = sum(auxArray)
# Instead present alleles, this shows missing ones
objetiveValue = sizeAlelleArray - x
exit = push!(exit, objetiveValue)
"""
Third objetive: maximize allele frequency
- sum of matrix B values
"""
y = 0.0
for i in 1:length(entry)
if entry[i] == 1
y = y + sum(data_matrixB[i,1:end])
end
end
exit = push!(exit,y)
"""
78
Fourth objetive: maximize heterozygozity -
matrix C values
Note: heterozygozity is a negative value -
greater heterozygozity is the lower
"""
objetiveValue = 0
for i in 1:length(entry)
if entry[i] == 1
objetiveValue = objetiveValue +
sum(data_matrixC[i, 1:end])
end
end
exit = push!(exit,objetiveValue)
"""
Fifth objetive: minimize HWE - matrix D values
"""
objetiveValue = 0
objetiveValue = abs(sum(data_matrixD.*entry))
exit = push!(exit,objetiveValue)
return exit
end
"""
Initializes a population, receives the population
size and returns an array of randomly created
Individuals
"""
function initPopulation(pop_size::Int)
population::Array{Individual} =[]
for i = 1:pop_size
79
population = push!(population, Individual(geneSize))
end
return population
end
########################################################
##
# NSGA-II in Julia
# Guilherme N. Ramos ([email protected])
# Gustavo Fernandes de Almeida ([email protected])
# Helper/utilitary functions.
##
"""
Computes the ranks and define fronts for individuals in
the population
"""
function fast_non_dominated_sort(population::Array{Individual})
fronts::Array{Front} = []
currentFront::Array{Individual} = []
nextFront::Array{Individual} = []
for p in population
p.S = []
p.n = 0
for q in population
if p != q
if dominates(p,q)
p.S = push!(p.S, q)
end
if dominates(q,p)
p.n = p.n + 1
end
80
end
end
if p.n == 0
p.rank = 1
currentFront = push!(currentFront, p)
end
end
fronts = push!(fronts, Front(currentFront))
i = 1 # Initialize front counter
while !isempty(fronts[i].individuals)
nextFront = []
for x in fronts[i].individuals
for y in x.S
y.n = y.n - 1
if y.n == 0
y.rank = i+1
nextFront = push!(nextFront,y)
end
end
end
i = i + 1
fronts = push!(fronts, Front(nextFront))
#fronts[i] = Q
end
return fronts
end
"""
Defines the crownding distance of a front’s individuals
"""
function crowding_distance_assigned(front::Array{Individual})
81
frontLen = length(front)
# Set every crowding distance to zero
for x in front
x.crowdingDistance = 0
end
n_objectives = length(front[1].fenotype)
for i in 1:n_objectives
sort!(front, lt = (x,y)-> x.fenotype[i] < y.fenotype[i],
rev = true)
front[1].crowdingDistance = front[end].crowdingDistance = Inf32
maximumValue = front[1].fenotype[i]
minimumValue = front[end].fenotype[i]
for y in 2:(frontLen - 1)
greaterNeighbor = front[y-1].fenotype[i]
lowerNeighbor = front[y+1].fenotype[i]
front[y].crowdingDistance = front[y].crowdingDistance +
(( greaterNeighbor - lowerNeighbor)/( maximumValue - minimumValue))
end
end
#println("--------------------------------------------------")
return
end
"""
Indicates if there should be a crossover.
"""
function should_crossover()
return CROSSOVER_PROBABILITY <= rand()
end
"""
Indicates if there should be a mutation.
82
"""
function should_mutate()
return MUTATION_PROBABILITY <= rand()
end
"""
Returns a random Individual from the given array.
"""
function random(population::Array{Individual})
index = rand(1:length(population))
return population[index]
end
"""
Randomly selects n Individuals from the given array
(with repetition).
"""
function random(n::Int, population::Array{Individual})
selected::Array{Individual} = []
for i in 1:n
push!(selected, random(population))
end
return selected
end
## ##
# Tournament #
## ##
"""
Tests if x and y are identical.
"""
function isequal(x::Individual, y::Individual)
83
return x.genotype == y.genotype
end
"""
Returns the winner of a tournament between x and y.
"""
function winner(x::Individual, y::Individual)
# Does x dominates y?
if (dominates(x,y))
return x
end
# Does y dominates x?
if (dominates(y,x))
return y
end
# No solution is definitely better
return random([x, y])
end
"""
Returns the winner of a binary tournament.
"""
function binary_tournament(contenders::Array{Individual})
x, y = contenders[1], contenders[2]
return isequal(x, y) ? random([x, y]) : winner(x, y)
end
"""
Tests if x dominates y.
In our scenario both fenotypes must be minimized
Example:
[2, 1] < [1,1]
[1, 2] < [1,1]
[1, 2] = [2,1] -> No answer is any better
84
"""
function dominates(x::Individual, y::Individual)
return x.fenotype[1] < y.fenotype[1] &&
x.fenotype[2] < y.fenotype[2] &&
x.fenotype[3] > y.fenotype[3] &&
x.fenotype[4] < y.fenotype[4] &&
x.fenotype[5] < y.fenotype[5]
end
"""
Doubles the size of the population.
"""
function expand_population(population::Array{Individual})
num_steps = length(population)/2
for i = 1:num_steps
first_parent = binary_tournament(random(2, population))
second_parent = binary_tournament(random(2, population))
# Clone
children = [Individual(parent.genotype)
for parent in [first_parent, second_parent]]
# Crossover
if should_crossover()
crossover(children[1], children[2])
#crossover(children)
end
# Mutation
for child in children
if should_mutate()
singleBitMutation(child)
end
85
end
append!(population, children)
end
end
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Variation functions
##
"""
Receives an array of integers and flips a bit
(0 to 1 or 1 to 0) in a random locus
"""
function singleBitMutation(gene::Array)
x = rand(1:length(gene))
if gene[x] > 0
gene[x] = 0
else
gene[x] = 1
end
end
"""
Receives a Individual and flips a bit in
a random locus
"""
function singleBitMutation(individualA::Individual)
x = rand(1:length(individualA.genotype))
if individualA.genotype[x] > 0
individualA.genotype[x] = 0
else
86
individualA.genotype[x] = 1
end
individualA.fenotype = initFenotype(individualA.genotype)
end
"""
Receives two vetors and switch their values from a
random locus (single point crossover)
"""
function crossover(geneA::Array{Int8} , geneB::Array{Int8})
if length(geneA) != length(geneB)
println("Erro during crossover")
return
end
i = rand(1:length(geneA))
v_aux = copy(geneA)
geneA[i:end] = geneB[i:end]
geneB[i:end] = v_aux[i:end]
end
"""
Crossover function version for Individual type
"""
function crossover(individualA::Individual,individualB::Individual)
if length(individualA.genotype) != length(individualB.genotype)
println("Erro during crossover")
return
end
i = rand(1:length(individualA.genotype))
v_aux = copy(individualA.genotype)
individualA.genotype[i:end] = individualB.genotype[i:end]
individualB.genotype[i:end] = v_aux[i:end]
individualA.fenotype = initFenotype(individualA.genotype)
individualB.fenotype = initFenotype(individualB.genotype)
87
end
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Display functions
##
"""
Prints a brief description of each Individual
in a population
"""
function printPopulation(p::Array{Individual})
for i = 1:length(p)
println("-----------------------------------")
println("Genotype:",repr(p[i].genotype))
println("Fenotype:",repr(p[i].fenotype))
println("S:",summary(p[i].S))
println("n:",repr(p[i].n))
println("Rank:", p[i].rank)
println("Crowding distance:",
p[i].crowdingDistance)
print("Press enter to continue...");
x = readline(STDIN)
end
end
"""
Prints in a file
"""
function createInterFile(p::Array{Individual})
outfile = open("inter.csv", "w")
88
for i = 1:length(p)
write(outfile, string(p[i].fenotype[1],",",
p[i].fenotype[2],",",
p[i].fenotype[3],",",
p[i].fenotype[4],",",
p[i].fenotype[5],"\n"))
end
close(outfile)
end
"""
Prints in a file and prompt some solutions, helps debugging
"""
function debug_printFile(p::Array{Individual})
outfile = open("exit.txt", "w")
for i = 1:length(p)
if p[i].fenotype[1] < 9 && p[i].fenotype[2] == 0
for j = 1:length(p[i].genotype)
if(p[i].genotype[j]>0)
write(outfile, string(j,"-") )
print(j, "-")
end
end
write(outfile,string("\n"))
print("\n")
end
end
close(outfile)
end
"""
Prints a summary of each element of a type
Taken from:
http://samuelcolvin.github.io/JuliaByExample/#Arrays
89
"""
function printsum(a)
# summary generates a summary of an object
println(summary(a), ": ", repr(a))
end
########################################################
##
# NSGA-II in Julia
# Gustavo Fernandes de Almeida ([email protected])
# Display functions - Plotting using Gadfly
##
using Gadfly
"""
Plots simple 2D graph that gives a summary of
a population distribution
"""
function plot_bee()
data = readcsv("inter.csv")
x1::Array{Int} = []
y1::Array{Int} = []
println("Plotting...")
for counter in 1:size(data)[1]
x1 = push!(x1, data[counter,1])
y1 = push!(y1, data[counter,2])
end
plot( x = x1,
y = y1,
Stat.xticks(ticks=[-1,0,1,2,3,4,5,6,7,
8,9,10,11,12,13,14,15,16,17,18,19,20,
21,22,23,24,25]),
Stat.yticks(ticks=[-5,0,5,10,15,20,25,
30,35,40,45,50,55]),
Guide.XLabel("How many areas"),
Guide.YLabel("How many missing alelles"),
90
Guide.title("NSGA-II’s results"),
Geom.beeswarm)
end
function plot_point()
data = readcsv("inter.csv")
x1::Array{Int} = []
y1::Array{Int} = []
println("Plotting...")
for counter in 1:size(data)[1]
x1 = push!(x1, data[counter,1])
y1 = push!(y1, data[counter,2])
end
plot( x = x1,
y = y1,
Stat.xticks(ticks=[-1,0,1,2,3,4,5,6,7,
8,9,10,11,12,13,14,15,16,17,18,19,20,
21,22,23,24,25]),
Stat.yticks(ticks=[-5,0,5,10,15,20,25,
30,35,40,45,50,55]),
Guide.XLabel("How many areas"),
Guide.YLabel("How many missing alelles"),
Guide.title("NSGA-II’s results"),
Geom.point)
end
91
Anexo II
Matriz A
Tabela II.1: Matriz A transposta.
Locus-Alelo CMT
ABMT
PCO
SMS
AMS
ATO
SMGO
LGO
ISP
MAMG
ENGO
STGO
AMG
PMG
PMS
PCMS
CMS
IGO
RAMT
RAGO
JGO
NTO
ART
OAQMS
CAMT
Bm164 - 178 2 0 0 2 6 2 1 0 0 1 4 4 0 1 1 0 8 0 1 1 0 0 0 9 2Bm164 - 176 0 0 0 19 2 0 0 0 12 25 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1Bm164 - 174 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2Bm164 - 170 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Bm164 - 168 0 0 1 0 12 0 0 0 0 0 0 1 0 0 0 1 6 0 0 0 0 0 0 0 0Bm164 - 165 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9Bm164 - 158 31 32 32 28 31 32 32 32 27 28 12 11 32 32 13 12 7 13 27 37 32 12 15 31 2Bm164 - 156 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 27DaE06 - 220 4 0 0 28 12 14 15 7 11 15 4 2 19 5 2 3 5 2 3 3 0 1 3 30 30DaE06 - 216 31 32 32 20 27 27 28 32 31 29 10 11 30 32 12 10 12 12 26 37 32 12 15 30 0DaE06 - 212 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0DaE12 - 222 3 2 0 0 1 8 8 0 14 12 0 0 0 0 1 1 0 2 8 12 9 0 1 2 0DaE12 - 220 32 32 32 31 31 32 30 32 21 30 12 12 32 32 10 12 13 13 27 37 31 12 15 31 28DaE12 - 219 0 16 0 0 2 0 0 0 6 1 0 0 0 0 1 4 7 8 0 1 12 0 0 0 0DaE12 - 218 0 0 0 3 8 0 0 0 0 0 0 0 0 0 9 3 3 0 0 0 0 0 0 0 12DaE12 - 216 0 0 0 0 3 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0DaE20 - 158 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3DaE20 - 156 8 5 20 5 2 20 3 21 0 24 0 2 4 0 5 2 2 7 17 13 16 4 3 19 4DaE20 - 154 32 29 28 30 27 23 32 29 31 32 12 12 32 32 12 10 13 13 20 36 31 12 15 25 24DaE20 - 146 0 11 0 1 11 0 0 0 0 0 0 0 0 0 0 4 2 0 2 1 0 0 0 2 4DaE34 - 124 7 12 2 0 5 1 0 0 0 0 0 0 0 0 0 0 0 5 13 0 4 0 0 0 0DaE34 - 122 5 26 30 12 9 2 2 18 16 18 4 9 22 2 0 2 3 8 0 18 23 1 6 0 0DaE34 - 120 0 0 0 0 0 0 0 0 0 1 2 1 1 1 4 0 11 0 1 3 0 0 0 1 0DaE34 - 118 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 2 0 0 0 0 0 0 0 0DaE34 - 116 11 8 11 5 3 6 3 20 15 0 2 1 0 0 1 0 0 0 7 6 0 7 6 21 0DaE34 - 114 16 0 0 5 1 15 15 3 2 11 1 3 0 0 3 1 0 0 5 6 7 1 5 2 23
Continua na próxima página
92
Tabela II.1 – continuando da última página
CMT
ABMT
PCO
SMS
AMS
ATO
SMGO
LGO
ISP
MAMG
ENGO
STGO
AMG
PMG
PMS
PCMS
CMS
IGO
RAMT
RAGO
JGO
NTO
ART
OAQMS
CAMT
DaE34 - 112 0 0 1 9 3 0 0 0 0 1 0 3 4 5 1 2 0 1 0 6 0 0 0 12 0DaE34 - 110 5 0 0 13 5 23 15 9 15 13 9 4 25 29 4 0 3 1 11 21 11 10 10 6 15DaE34 - 108 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 1 0DaE34 - 106 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0DaE34 - 104 0 0 0 0 11 0 12 4 0 0 0 0 0 1 4 4 0 0 0 0 0 0 0 1 0DaE41 - 150 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 1 0 0 0 0 0 0 0DaE41 - 148 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 2 0 0 0 1DaE41 - 146 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0DaE41 - 142 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7DaE41 - 138 0 0 0 4 8 0 0 0 0 0 0 0 0 0 3 0 0 0 0 8 0 0 0 0 2DaE41 - 136 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 4DaE41 - 134 0 0 0 4 1 4 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0DaE41 - 132 0 0 1 2 3 15 0 0 0 4 0 4 1 0 0 0 0 0 0 0 0 0 0 0 5DaE41 - 130 5 18 11 0 1 0 28 24 0 5 7 1 0 5 6 0 2 10 15 7 20 9 3 8 0DaE41 - 128 0 1 0 10 7 2 0 0 0 0 0 0 0 0 0 1 1 0 1 10 0 0 0 0 9DaE41 - 126 31 15 21 0 0 4 16 3 2 19 6 10 32 0 0 0 0 0 0 0 1 0 1 10 0DaE41 - 124 0 0 0 19 3 11 5 18 30 0 0 0 0 31 0 1 3 0 15 0 20 10 12 22 10DaE41 - 122 0 0 0 8 16 3 0 0 1 12 1 0 0 0 1 12 9 2 5 19 0 0 0 1 0DaE41 - 120 0 0 5 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 4 5 3 0 0 1 5DaE41 - 118 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 2DaE63 - 214 0 0 16 10 4 0 0 5 6 3 0 0 0 0 0 1 0 2 0 0 22 0 0 17 15DaE63 - 210 21 7 17 9 2 28 4 23 19 16 10 8 12 15 1 1 2 11 10 26 13 12 10 1 0DaE63 - 208 26 31 19 25 32 16 32 12 13 28 8 11 28 25 13 12 13 6 25 36 24 6 12 26 22DaE67 - 176 32 30 32 31 32 22 32 32 31 32 12 12 14 6 11 13 13 13 27 37 32 12 15 31 30DaE67 - 170 5 8 0 6 5 27 15 0 0 0 5 1 26 32 7 2 0 0 0 0 0 0 1 2 3DaE46 - 253 28 12 16 2 2 22 28 12 14 27 10 7 1 6 3 10 4 4 6 11 25 12 15 10 5DaE46 - 250 4 22 17 29 30 10 0 20 17 5 2 5 29 25 9 2 9 9 20 26 7 0 0 21 25DaE46 - 247 8 19 2 13 23 22 11 6 13 6 8 6 4 8 2 3 13 8 6 32 27 7 13 16 5DaE46 - 244 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
93
Anexo III
Matriz B
94
Tabela III.1: Matriz B transposta. Valores multiplicados por 104 para melhor visualização.
Locus-Alelo CMT
ABMT
PCO
SMS
AMS
ATO
SMGO
LGO
ISP
MAMG
ENGO
STGO
AMG
PMG
PMS
PCMS
CMS
IGO
RAMT
RAGO
JGO
NTO
ART
O
AQMS
CAMT
Bm164-178 1.8 0 0 1.7 4.8 1.6 0.8 0 0 0.8 23.6 23.6 0 1 5.3 0 36.6 0 1.2 0.6 0 0 0 7.5 2Bm164-176 0 0 0 15.9 1.6 0 0 0 11.2 19.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1Bm164-174 0 0.8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2Bm164-170 4.4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Bm164-168 0 0 0.9 0 9.6 0 0 0 0 0 0 5.9 0 0 0 5.7 27.5 0 0 0 0 0 0 0 0Bm164-165 0 0 0 0.8 1.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8.9Bm164-158 27.5 27.1 28.9 23.4 24.9 25.6 26.4 27.6 25.1 22 70.9 65 28.7 30.8 68.5 68.9 32.1 65.8 32.9 21.7 24.4 70.9 56.8 25.7 2Bm164-156 0 0 0 1.7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 26.8DaE06-220 3.6 0 0 23.4 9.6 11.2 12.4 6 10.2 11.8 23.6 11.8 17.1 4.8 10.5 17.2 22.9 10.1 3.7 1.8 0 5.9 11.4 24.9 29.8DaE06-216 27.5 27.1 28.9 16.7 21.7 21.6 23.1 27.6 28.8 22.8 59.1 65 26.9 30.8 63.2 57.4 54.9 60.7 31.7 21.7 24.4 70.9 56.8 24.9 0DaE06-212 0 0 0 0 0.8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0DaE12-222 2.7 1.7 0 0 0.8 6.4 6.6 0 13 9.4 0 0 0 0 5.3 5.7 0 10.1 9.7 7 6.9 0 3.8 1.7 0DaE12-220 28.4 27.1 28.9 25.9 24.9 25.6 24.7 27.6 19.5 23.6 70.9 70.9 28.7 30.8 52.7 68.9 59.5 65.8 32.9 21.7 23.6 70.9 56.8 25.7 27.8DaE12-219 0 13.6 0 0 1.6 0 0 0 5.6 0.8 0 0 0 0 5.3 23 32.1 40.5 0 0.6 9.1 0 0 0 0DaE12-218 0 0 0 2.5 6.4 0 0 0 0 0 0 0 0 0 47.4 17.2 13.7 0 0 0 0 0 0 0 11.9DaE12-216 0 0 0 0 2.4 0 0 0 0 0 0 0 0 0 5.3 0 0 0 0 0 0 0 0 0 0DaE20-158 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3DaE20-156 7.1 4.2 18.1 4.2 1.6 16 2.5 18.1 0 18.8 0 11.8 3.6 0 26.3 11.5 9.2 35.4 20.7 7.6 12.2 23.6 11.4 15.8 4DaE20-154 28.4 24.6 25.3 25.1 21.7 18.4 26.4 25 28.8 25.1 70.9 70.9 28.7 30.8 63.2 57.4 59.5 65.8 24.4 21.1 23.6 70.9 56.8 20.7 23.8DaE20-146 0 9.3 0 0.8 8.8 0 0 0 0 0 0 0 0 0 0 23 9.2 0 2.4 0.6 0 0 0 1.7 4DaE34-124 6.2 10.2 1.8 0 4 0.8 0 0 0 0 0 0 0 0 0 0 0 25.3 15.8 0 3 0 0 0 0DaE34-122 4.4 22 27.1 10 7.2 1.6 1.6 15.5 14.9 14.1 23.6 53.2 19.8 1.9 0 11.5 13.7 40.5 0 10.6 17.5 5.9 22.7 0 0DaE34-120 0 0 0 0 0 0 0 0 0 0.8 11.8 5.9 0.9 1 21.1 0 50.4 0 1.2 1.8 0 0 0 0.8 0DaE34-118 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 28.7 9.2 0 0 0 0 0 0 0 0DaE34-116 9.8 6.8 9.9 4.2 2.4 4.8 2.5 17.3 13.9 0 11.8 5.9 0 0 5.3 0 0 0 8.5 3.5 0 41.4 22.7 17.4 0DaE34-114 14.2 0 0 4.2 0.8 12 12.4 2.6 1.9 8.6 5.9 17.7 0 0 15.8 5.7 0 0 6.1 3.5 5.3 5.9 18.9 1.7 22.8DaE34-112 0 0 0.9 7.5 2.4 0 0 0 0 0.8 0 17.7 3.6 4.8 5.3 11.5 0 5.1 0 3.5 0 0 0 10 0DaE34-110 4.4 0 0 10.9 4 18.4 12.4 7.8 13.9 10.2 53.2 23.6 22.4 27.9 21.1 0 13.7 5.1 13.4 12.3 8.4 59.1 37.9 5 14.9DaE34-108 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4.6 0 0 0.8 0DaE34-106 0 0 0 0 0 0 9.9 0 0 0 0 0 0 0 0 0 9.2 0 0 0 0 0 0 0 0DaE34-104 0 0 0 0 8.8 0 9.9 3.5 0 0 0 0 0 1 21.1 23 0 0 0 0 0 0 0 0.8 0DaE41-150 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15.8 0 0 5.1 0 0 0 0 0 0 0DaE41-148 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3.5 1.5 0 0 0 1DaE41-146 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.2 0 0 0 0 0 0
Continua na próxima página
95
Tabela III.1 – continuando da última páginaCMT
ABMT
PCO
SMS
AMS
ATO
SMGO
LGO
ISP
MAMG
ENGO
STGO
AMG
PMG
PMS
PCMS
CMS
IGO
RAMT
RAGO
JGO
NTO
ART
O
AQMS
CAMT
DaE41-142 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.9DaE41-138 0 0 0 3.3 6.4 0 0 0 0 0 0 0 0 0 15.8 0 0 0 0 4.7 0 0 0 0 2DaE41-136 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5.3 0 0 0 0 0 0 0 0 0 4DaE41-134 0 0 0 3.3 0.8 3.2 0 0 0 0 0 0 0 0 0 0 0 5.1 0 0 0 0 0 0 0DaE41-132 0 0 0.9 1.7 2.4 12 0 0 0 3.1 0 23.6 0.9 0 0 0 0 0 0 0 0 0 0 0 5DaE41-130 4.4 15.2 9.9 0 0.8 0 23.1 20.7 0 3.9 41.4 5.9 0 4.8 31.6 0 9.2 50.6 18.3 4.1 15.2 53.2 11.4 6.6 0DaE41-128 0 0.8 0 8.4 5.6 1.6 0 0 0 0 0 0 0 0 0 5.7 4.6 0 1.2 5.9 0 0 0 0 8.9DaE41-126 27.5 12.7 19 0 0 3.2 13.2 2.6 1.9 14.9 35.5 59.1 28.7 0 0 0 0 0 0 0 0.8 0 3.8 8.3 0DaE41-124 0 0 0 15.9 2.4 8.8 4.1 15.5 27.9 0 0 0 0 29.8 0 5.7 13.7 0 18.3 0 15.2 59.1 45.5 18.2 9.9DaE41-122 0 0 0 6.7 12.9 2.4 0 0 0.9 9.4 5.9 0 0 0 5.3 68.9 41.2 10.1 6.1 11.1 0 0 0 0.8 0DaE41-120 0 0 4.5 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 4.9 2.9 2.3 0 0 0.8 5DaE41-118 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10.5 0 0 0 0 0 0 5.9 0 0 2DaE63-214 0 0 14.5 8.4 3.2 0 0 4.3 5.6 2.4 0 0 0 0 0 5.7 0 10.1 0 0 16.8 0 0 14.1 14.9DaE63-210 18.6 5.9 15.4 7.5 1.6 22.4 3.3 19.9 17.7 12.6 59.1 47.3 10.8 14.4 5.3 5.7 9.2 55.7 12.2 15.2 9.9 70.9 37.9 0.8 0DaE63-208 23.1 26.3 17.2 20.9 25.7 12.8 26.4 10.4 12.1 22 47.3 65 25.1 24 68.5 68.9 59.5 30.4 30.5 21.1 18.3 35.5 45.5 21.6 21.8DaE67-176 28.4 25.4 28.9 25.9 25.7 17.6 26.4 27.6 28.8 25.1 70.9 70.9 12.6 5.8 58 74.6 59.5 65.8 32.9 21.7 24.4 70.9 56.8 25.7 29.8DaE67-170 4.4 6.8 0 5 4 21.6 12.4 0 0 0 29.6 5.9 23.3 30.8 36.9 11.5 0 0 0 0 0 0 3.8 1.7 3DaE46-253 24.9 10.2 14.5 1.7 1.6 17.6 23.1 10.4 13 21.2 59.1 41.4 0.9 5.8 15.8 57.4 18.3 20.2 7.3 6.4 19.1 70.9 56.8 8.3 5DaE46-250 3.6 18.6 15.4 24.2 24.1 8 0 17.3 15.8 3.9 11.8 29.6 26 24 47.4 11.5 41.2 45.5 24.4 15.2 5.3 0 0 17.4 24.8DaE46-247 7.1 16.1 1.8 10.9 18.5 17.6 9.1 5.2 12.1 4.7 47.3 35.5 3.6 7.7 10.5 17.2 59.5 40.5 7.3 18.8 20.6 41.4 49.2 13.3 5DaE46-244 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.2 0 0 0 0 0 0
96
