Shell x Listas x Monte Carlo

Muitas vezes é interessante guardas listas de sítios e sortear somente dentro desta lista para evitar gerar números aleatórios que não serão usados. É fácil usar listas em uma dimensão: quando um sítio sai da lista ele pode ser substituído pelo ultimo da fila e esta perde um elemento.

Tudo muito bom, tudo muito legal, mas há um preço a pagar se estamos interessados no comportamento dinâmico... um não, dois !

Gerando números aleatórios com a biblioteca VSL (MKL)

Exemplo simples do uso dos geradores de números aleatórios da biblioteca MKL,

#include <stdio.h>
/////////////////////////////// MKL GSL Mersenne Twister
#include <mkl_vsl.h>

#define BRNG    VSL_BRNG_MT19937 
#define METHOD  0
////////////////////////////////

main ()
{

  int N = 4, sqrtN = 2 , i , j;
  int seed;
  seed = 345;
  
  float ranfloat[sqrtN]; //cada entrada corresponde a um numero aleatorio
  VSLStreamStatePtr stream; 
  vslNewStream (&stream, BRNG, seed); //semeando o gerador
  
  for (i = 0; i < sqrtN; i++)
    {
    vsRngUniform (METHOD, stream, sqrtN, ranfloat, 0.0, 1.0); //carregando sqrtN numeros
    for (j = 0; j < sqrtN; j++)
      {
      printf ("%f\n", ranfloat[sqrtN]);
      }
    }

  return (0);

}

Para compilar basta linkar as bibliotecas adequadas:

ia32,  icc exemplo.c  -lmkl_intel -lmkl_sequential -lmkl_core  -o exemplo.bin

intel64, icc exemplo.c -lmkl_intel_lp64 -lmkl_sequential -lmkl_core  -o exemplo.bin

Pontos importante:
os números são obtidos em stream, a idéia é carrega um vetor a cada chamada da função vsRngUniform, no exemplo carregamos sqrtN números entre 0 e 1 (float) distribuídos uniformente. O desempenho é melhor para sqrtN > 1000.

Outras funções: 
vdRngUniform  -  double
viRngUniform  - int

Outros geradores, funções e mais detalhes podem ser obtidos na documentação da biblioteca Vector Statistical Library VSL

Geradores de números pseudo-aleatórios

Geradores de números pseudo-aleatórios GNPA são uma caixa preta, em geral usamos as indicações dos estudantes mais antigos e dos professores. Provavelmente, graças a esse mecanismo, os GNPA mais usados são o ran2 ou o rand, e a implementação destes é encontrada no clássico Numerical Recipes.

O cerne das simulações de Monte-Carlo são justamente os GNPA e dificilmente nossos programas passarão fazendo mais tempo algo que não seja gerar números aleatórios. Por isso é importante conhecer e se atualizar quanto este assunto.

Embora seja importante entender um pouco da minúcias dos GNPA eu não vou nem tocar nesse assunto, para mim eles permaneceram uma caixa preta por muito tempo. Bom..., as atualidades ? Ah, essas são bem mais fáceis. 

O GNPA padrão nos dias de hoje o Mersenne-Twister, introduzido em 1997 http://en.wikipedia.org/wiki/Mersenne_twister (alguem tem que escrever o artigo na wikipedia portuguesa, nem que seja tradução). Este é provavelmente o algoritmo mais popular, ele possui um período de 2^19937 − 1, passou em vários testes e é superior aos ran2,3,4 .

Outro ponto de fundamental importância é a sua performance.

Como implementar o algoritmo:

Eu sugiro fortemente a utilização da biblioteca GSL ( Gnu Scientific Library ) http://www.gnu.org/software/gsl/ e o manual http://www.gnu.org/software/gsl/manual , essa biblioteca já esta instalada no abax e você pode encontra-la nos repositórias da sua distribuição Linux.

E muito fácil usar, segue um exemplo

#include < stdio.h >
#include < gsl/gsl_rng.h >
int main ()
  {
  int i;
  unsigned long a;
/////////////////////////////////////////////////////
  const gsl_rng *r;
  r = gsl_rng_alloc(gsl_rng_mt19937);
////////////////////////////////////////////////////
  gsl_rng_set(r, 123); // semente 123
   for (i = 0; i < 100000000; i++)
    {
    printf ("%u\n",  gsl_rng_get(r));
    }
  printf ("\n");
  printf ("%lu %lu\n", gsl_rng_max(r), gsl_rng_min(r) );
  return 0;
}

Compilar com: gcc rng.c -lgsl -lgslcblas -lm

Comentários:

"const gsl_rng *r;
r = gsl_rng_alloc(gsl_rng_mt19937); "
- Define o nome do gerador, no caso é "r",  e o GNPA usado ( mt19937 = Mersenne Twisted ).

"gsl_rng_set (r, 123);" - Semente (123) do gerador r.

"gsl_rng_get (r);" - Retorna um número inteiro (unsigned long) pseudo-aleatório

"gsl_rng_max(r) e gsl_rng_min(r)" - Valor máximo e mínimo do gerador.


Outras dicas:

"gsl_rng_uniform(r);" - Retorna um real (double) pseudo-aleatório o intervalo [0,1)
"gsl_rng_uniform_int(r,j);" Retorna um número inteiro (unsigned long) pseudo-aleatório entre [0, j-1]

Existe uma lista de extensa de implementações de GNPA nesta biblioteca, inclusive o ran2, basta substituir o "mt19937" pelo algoritmo desejado.

Veja a lista de performance ( http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generator-Performance.html ):

     1754 k ints/sec,    870 k doubles/sec, taus
     1613 k ints/sec,    855 k doubles/sec, gfsr4
     1370 k ints/sec,    769 k doubles/sec, mt19937
      565 k ints/sec,    571 k doubles/sec, ranlxs0
      400 k ints/sec,    405 k doubles/sec, ranlxs1
      490 k ints/sec,    389 k doubles/sec, mrg
      407 k ints/sec,    297 k doubles/sec, ranlux
      243 k ints/sec,    254 k doubles/sec, ranlxd1
      251 k ints/sec,    253 k doubles/sec, ranlxs2
      238 k ints/sec,    215 k doubles/sec, cmrg
      247 k ints/sec,    198 k doubles/sec, ranlux389
      141 k ints/sec,    140 k doubles/sec, ranlxd2
    
     1852 k ints/sec,    935 k doubles/sec, ran3
      813 k ints/sec,    575 k doubles/sec, ran0
      787 k ints/sec,    476 k doubles/sec, ran1
      379 k ints/sec,    292 k doubles/sec, ran2


Uma pequena explicação sobre as características de cada algoritmo é encontrada no manual http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html


Using GSL to calculate Nonsymmetric matrices eigenvalues

Here is a example of using GSL - Gnu Scientific Library - to calculate the eigenvalues of a nonsymmetric matrix.
This is a simple example write in C. GSL uses BLAS and can be optimized with ATLAS, so can be used as a real tool for heavy calculations. To get know more about GSL visit http://www.gnu.org/software/gsl/ or the manual page.

Hint: Beware when using GSL or our own linear algebra routine, remember that greats matrices requires greats responsibilities, so keep track of precision errors. Find precision errors will be a great excuse to know more about multiple precision GMP .

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
int main (void)
        {
        int a=2; //matrix a x a
        int i,j;
        double cg[a][a];

        // Defining matrix
        cg[0][0] = 2.;
        cg[0][1] = 3.;
        cg[1][0] = -1.;
        cg[1][1] = 1.;

        gsl_matrix *TOTAL = gsl_matrix_alloc (a,a); //Alloc GSL matrix TOTAL

        //Setting TOTAL
        for (i = 0; i < a; i++)
                {
                for (j = 0; j < a; j++)
                        {
                        gsl_matrix_set (TOTAL, i, j, cg[i][j]);
                        }
                }


        gsl_vector_complex *eval = gsl_vector_complex_alloc (a); //Vector whose entries are the eigenvalues
        gsl_eigen_nonsymm_workspace *K = gsl_eigen_nonsymm_alloc (a); //alloc work space for GSL nonsymmetric calculations
        gsl_eigen_nonsymm (TOTAL,eval,K);
        gsl_eigen_nonsymm_free (K);
        gsl_matrix_free (TOTAL);


        for (i = 0; i < a; i++)  /* OUT OF RANGE ERROR */
                {
                gsl_complex eval_i = gsl_vector_complex_get (eval, i);
                printf ("%g + %g*i\n",GSL_REAL(eval_i),GSL_IMAG(eval_i) );
                }
return 0;
}