q-Voter em networks
Junho-julho/2026, Mestrado, TF.
Implementação do modelo em redes com distribuição de graus \(p(k)\) (Poisson, lei de potência, \(\dots\)).
Mean-field e q-Voter
O ponto de partida aqui, então, é o artigo:
📄 Mean-Field Analysis of the q-Voter Model on Networks
Paolo Moretti, Suyu Liu, Claudio Castellano, and Romualdo Pastor-Satorras
Journal of Statistical Physics 151, 113–130 (2013).
Neste artigo, os autores implementam a teoria de campo médio num grafo completo e com a teoria heterogênea, HMF. Nos dois casos, a teoria prevê duas fases, para (fase de fragmentação da opinião média) e ferromagnética (consenso da opinião média), bem como um região de biestabilidade, que, dependendo da condição inicial, leva ao consenso ou à fragmentação.
A curva que define a fronteira crítica entre a fase paramagnética e ferromagnética é dada por \[ \overline{\varepsilon_c}(\phi) = \frac{\phi^2+3}{2\phi^2+14}\,, \tag{1}\] onde \(\phi\) é a magnetização (opinião) média da rede. Para a teoria HMF, sendo \(\langle k\rangle\) o grau médio da rede, a fronteira é dada por \[ \overline{\varepsilon_c}(\phi) = \frac{1}{2} - \frac{2\langle k\rangle(\langle k\rangle-2)}{\phi^2(\langle k\rangle-2) \tag{2}\]
Algumas conclusões centrais deste artigo: - As teorias de campo médio prevem biestabilidade para o q-voter para \(q\geq4\), cuja fronteira é dada pelas equações {ref}fronteira_critica e {ref}fronteira_critica; - A taxa do modelo em redes heterogêneas fica definida para um nó de grau \(k\) com \(n\) vizinhos em um estado diferente por: \[f_k(n;q,\varepsilon) = \left(\frac{n}{k}\right)^q + \varepsilon\left[1 - \left(\frac{n}{k}\right)^q - \left(\frac{k-n}{k}\right)^q\right]\,;\] - Em RRNs com \(K=4\), a biestabilidade não é observada, como previsto pelas teorias de campo médio, por conta da esparsidade da rede; - Questões se RRNs muito densas e redes heterogêneas, como as em lei de potência, recuperam a biestabilidade permanecem abertas no artigo.
Nos interessa, aqui, o caso \(q=4\), o menor valor de \(q\) em que a biestabilidade é observada na teoria de campo médio. Procuramos fazer um diagrama de fases semelhante ao da Figura {ref}biestabilidade, onde, para um dado valor \(\varepsilon\), examinamos a dinâmica varrendo as condições iniciais da magnetização, \(\phi_0\in[0,1]\).
Implementação computacional: benchmarks e resultados
Nas subseções que seguem, explico detalhadamente de que forma as redes complexas e a dinâmica (em tempo contínuo) foram implementadas. Também apresento benchmarks sobre a coerência e exatidão dos algoritmos, bem como as otimizações (lista dinâmica e paralelismo).
### Redes com diferentes distribuições de graus
A implementação da rede utiliza de tipos derivados em Fortran. Em suma, uma rede corresponde a um objeto com algumas propriedades intrínsecas, a saber: - adjacency_list (variável alocável): guarda o grau, \(k_i\), e a lista de vizinhos de cada nó \(i\) da rede (coleção de \(N\) vetores alocáveis, cada um de tamanho \(k_i\)); - gen_network (sub-rotina): gera a lista de graus da rede a partir de uma dada distribuição (power law, Poisson ou constante); - build_network (sub-rotina): aplicação do modelo configuracional a partir de uma dada distribuição de graus e o tamanho da rede N; - report (função): verifica os possíveis erros e inconsistências que podem aparecer com o modelo configuracional.
A grande vantagem desta abordagem é a capacidade de gerar redes totalmente independentes sem acúmulo de variáveis (ou uso de dummies), compacticidade do código e possibilidade de paralelização com OpenMP (basta passar a rede, e todas suas propriedades, como variável privada de uma thread), sem custo computacional extra.
Para gerar a distribuição de graus, dois algoritmos foram utilizados. No caso de uma distribuição power-law, o método mais interessante é a inversão de CDF com aproximação para dados discretos (acrescentar citação do Clauset); isto é, considere um número aleatório \(u^{\text{r}}\in[0,1)\), então o número inteiro aleatório \(x^{\text{r}}\), dado por: \[ x^{\text{r}} = \left\lfloor \big(x_{\min}-0.5\big)\big(1-u^{\text{r}}\big)^{\frac{1}{\gamma-1}}\right\rfloor\,, \tag{3}\] é um número distribuído em lei de potência \(p(x;\gamma)\) no intervalo \([x_{\min},\infty)\). Para controlar o surgimento de outliers e/ou o grau máximo da rede em lei de potência, pode-se atribuir um determinado cut-off: estrutural \(\max(k) = \sqrt{N}\), natural \(\max(k) = N^{\frac{1}{\gamma-1}}\) ou rígido \(\max(k) = N^{\frac{1}{\gamma}}\)
Para a de Poisson, o algoritmo de D. Knuth (ou sua variante, Junhao, acrescentar citações) é ideal. Ele parte da ideia que um número de Poisson conta o número de eventos que ocorrem numa janela temporal, em que o intervalo de tempo entre eventos sucessivos é dado por uma distribuição exponencial. Em outras palavras, é gerar números exponencialmente distribuídos com média \(1/\lambda\) de tal forma que o número \(N^{\text{r}}=m-1\), em que \(m\) é o menor número tal que a sequência gerada é maior ou igual a um (isto é, \(x_1^{\text{r}}+x_2^{\text{r}}+\cdots+x_m^{\text{r}}\geq 1\)), é um número com distribuição de Poisson com média \(\lambda\). Por razões de consistência, a distribuição gerada é truncada à esquerda, rejeitando graus inferiores a um dado \(k_{\min}\).
O modelo configuracional segue o algoritmo abaixo: - Dada a sequência de graus gerada, cria-se um vetor de “pontas soltas” (stubs) de tamanho \(\sum k_i\) em que cada nó \(i\) aparece (sequencialmente) \(k_i\) vezes; - Para evitar correlações na geração sequencial (seeds colision), o vetor de stubs é embaralhado com o Algoritmo P de Fisher-Yates (citar o livro do Knuth); - O vetor é percorrido emparelhando stubs adjacentes (índices \(j\) e \(j+1\)) para formar uma possível nova aresta entre os nós: - Caso a possível conexão for um auto-laço ou já existir, ela é rejeitada e o algoritmo continua: - Em caso de rejeição, um stub é sorteado (uniforme) dentre aqueles ainda não processados; - Caso falhe novamente, uma busca linear procura pelo primeiro stub válido disponível. Se nenhuma reatribuição for possível, o par é descartado; - Se a conexão for válida, a lista de adjacência de ambos os nós é atualizada, contabilizando a aresta final, e o processo avança até o esgotamento do vetor de stubs.
O método de sorteio e busca linear (two-stage cascade) serve para tentar ao máximo (para uma dada seed) encontrar conexões válidas para aquela distribuição. A {ref}benchmark_dist_efetiva_redes mostra que a distribuição efetiva de graus (a distribuição de graus real da rede, após a aplicação do modelo configuracional) de fato segue as distribuições usadas para a amostragem dos graus, com boa convergência mesmo para redes com \(N=10^4\) nós.
:name: benchmark_dist_efetiva_redes
:width: 100%
:align: center
Distribuições *efetivas* (reais) das redes geradas para dois tamanhos, $N=5\cdot10^5$ e $N=10^4$, em todo caso com $k_{\min}=4$. As linhas sólidas representam as distribuições teóricas com estes parâmetros, mostrando a concordância da distribuição efetiva obtida com o modelo configuracional e a "distribuição de entrada".Finalmente, para verificar a autoconsistência das redes sintéticas e medir os desvios da lista de graus induzidos pelo método, a função de report faz uma grande checagem estrutural e estatística da topologia gerada. A verificação é dividida em duas etapas:
- Autoconsistência Topológica:
- Teorema do Grau: o número de arestas tem que ser igual à metade da soma dos graus efetivos dos nós;
- Nós isolados: Contabiliza a existência de sítios com grau nulo (\(k=0\));
- Fragmentação da rede: Utilizando o algoritmo BFS, o código determina o tamanho da componente gigante (CG), de modo que, caso ela não equivalha à totalidade da da rede (sinalizando fragmentação), um alerta é emitido para que a componente principal seja isolada/extraída antes de iniciar a dinâmica.
- Consistência Estatística (Teoria vs. Efetiva):
- Arestas rejeitadas: Compara o número de arestas originais sorteadas com o número de arestas efetivamente conectadas, indicando a taxa de descarte exigida para manter a rede puramente simples;
- Momentos da distribuição: Calcula o grau máximo (\(k_{\max}\)), o grau médio (\(\langle k \rangle\)) e o segundo momento matemático (\(\langle k^2 \rangle\)) da lista original e compara com a rede conectada. Uma tolerância estatística de \(1/\sqrt{N}\) é adotada para classificar a confiabilidade da topologia gerada.
Estes relatórios on-the-fly servem para garantir robustez e validade ao algoritmo de redes, evitando que desvios induzidos pelo método enviesem os resultados.
Para várias redes geradas, nenhum desvio significativo foi encontrado. Assim, a robustez e validade deste algoritmo ficam comprovadas. A tabela abaixo reúne estes resultados comparativos para os casos de interesse.
| Distribuição | Parâmetros | Desvios da rede gerada |
|---|---|---|
| Power-law (Scale-Free) | \(\gamma \in \{2.3, 2.7, 3.5\}\), \(k_{\min}=4\) | “Insignificantes” Frequentemente 0 (ou \(\sim 10^{-5}\%\)) de descarte; CG abrange todos os nós; \(\langle k \rangle\), \(\langle k^2 \rangle\) e \(k_{\max}\) efetivos são idênticos aos da lista amostrada. |
| Poisson (Erdős-Rényi) | \(\langle k \rangle_{\text{teor}} \in \{6.0, 8.5, 10.0\}\), \(k_{\min}=4\) | Nulos. 0 arestas descartadas; CG abrange todos os nós; momentos estatísticos perfeitamente conservados. |
| Regular (RRN) | \(K \in \{4, 6, 10\}\) | Nulos. Conservação exata da topologia regular; 0 arestas descartadas; CG engloba toda a rede. |
### Dinâmica em tempo contínuo
Para utilizar o algoritmo de Gillespie, o método de aceitação-rejeição (eventos reais e nulos) com lista dinâmica é a forma otimizada. A referência principal para implementação do algoritmo foi: > 📄 Dynamic sampling from a discrete probability distribution with a known distribution of rates
> Federico D’Ambrosio, Hans L. Bodlaender, and Gerard T. Barkema
> Computational Statistics 37, 1203–1228 (2022).
A ideia é considerar apenas os nós que podem flipar. Basta definir uma lista que contém todos os \(N_{\text{ativos}}\) nós ativos da rede (isto é, aqueles com \(f_k>0\)) e atualizá-la dinamicamente durante a evolução. Em suma: - A um dado tempo \(t\), sorteia-se aleatoriamente um nó \(i\) exclusivamente dentre os \(N_\text{ativos}\) sítios com taxa não-nula; - Aceitação-rejeição padrão: dado um número aleatório uniforme \(r^\text{r}\in[0,1)\), o spin flipa caso \(f_{k_i}(n_i) \geq r^\text{r} \cdot \max_f\); - Caso o estado seja alterado, a vizinhança local é atualizada. Sítios vizinhos que passam a ter \(0\) vizinhos discordantes são removidos da lista dinâmica (tornam-se inativos), enquanto sítios que adquirem discordância são adicionados à lista; - O tempo físico do sistema avança independentemente do flip ter sido aceito ou rejeitado (evento nulo), considerando \(u^\text{r}\in[0,1)\), por: \[t \leftarrow t - \frac{\ln(1-u^\text{r})}{N_\text{ativos} \cdot \max_f}\,.\]
No pior dos casos, a lista dinâmica tem complexidade temporal igual ao método ingênuo de aceitação-rejeição (onde o sorteio do nó \(i\) é feito sobre toda a rede). Mas em dinâmicas próximas ao estado de consenso, o tamanho da lista decresce rapidamente, de forma que o uso da lista dinâmica reduz a complexidade temporal. Outros métodos de otimização, a princípio, não geram melhorias no tempo computacional.
:name: benchmark_equiv
:width: 80%
:align: center
Equivalência entre o algoritmo de Gillespie otimizado com a lista dinâmica e o algoritmo em passos discretos de Markov.A Figura {ref}benchmark_equiv mostra a equivalência (estatística) das implementações em tempo discreto e tempo contínuo com lista dinâmica. As Figuras abaixo mostram o ganho computacional com o método em tempo contínuo.
 |
 |
### Amostragem dos resultados
Para obtenção dos resultados, a condição de parada da simulação é dada por um tempo máximo numericamente igual ao tamanho da rede Gillespie_time = t_max=real(N) ou à obtenção do consenso, N_ativos = 0. Isto é feito para diferenciar claramente as fases de fragmentação e consenso. Na fase paramagnética o consenso é obtido em tempos muito grandes (\(\tau\sim\exp(N)\)), de modo que este “tempo máximo” serve para relevar o efeito de tamanho finito. De qualquer forma, a análise de tamanho finito ainda não foi realizada.
Para reproduzir a {ref}refbiestabilidade, foram 50 amostras (feitas com uma única rede) de cada ponto no espaço \((\phi_0,\varepsilon)\). Aqui, \(\phi_0\) entra como condição inicial: \((\phi_0+1)/2\) é a proporção (probabilidade) de agentes no estado \(+1\). Os diagramas abrangem 123 pontos distribuídos linearmente para \(\varepsilon\in[0,0,0.3]\) e 101 para \(\phi_0\in[-0.95,0.95]\) (totalizando 12423 pontos * 50 amostras = 621150 runs para cada rede).
Para adiantar a obtenção de amostras, a ferramenta OpenMP foi utilizada, utilizando somente núcleos de desempenho das máquinas do GISC. Para realizar o teste de benchmarking, 100 amostras (50 na fase paramagnética e 50 na fase ferromagnética) foram rodadas sequencialmente (em único núcleo) e com/sem hyper-threading na máquina Weiss. Como os processadores das máquinas Ising, Parisi, Fisher e Weiss seguem a mesma topologia (verificado com lstopo), o resultado percentual da tabela abaixo se mantém. O atraso do hyper-threading se explica porque a memória que a dinâmica ocupa é maior que a metade da memória disponível no cache L2: matriz de estados: 80 KB, matriz de taxas: 80 KB, lista de adjacência: 781 KB e lista dinâmica: 160KB, memória total de aproximadamente 1,1 MB, o que inviabiliza alocar duas threads simultâneas disputando o mesmo cache L2 de 2 MB sem causar cache misses.
| Execução | Nº de Threads | Speedup (Ganho) | Eficiência Paralela | Tempo de Execução (s) |
|---|---|---|---|---|
| Sequencial Puro | 1 | 1,00x | 100 % | 351.013 |
| Sem Hyper-Threading | 8 | 6,05x | 75,7 % | 57.954 |
| Com Hyper-Threading | 16 | 5,74x | 35,9 % | 61.087 |
Resultados preliminares
:name: label
:width: 100%
:align: center
legenda:name: label
:width: 100%
:align: center
legenda:name: label
:width: 50%
:align: center
legenda:name: label
:width: 100%
:align: center
legenda:name: label
:width: 100%
:align: center
legenda:name: label
:width: 100%
:align: center
legenda