Quão lento é realmente o Python? (Ou quão rápido é o seu idioma?)

149

Eu tenho esse código que eu escrevi em Python / NumPy

from __future__ import division
import numpy as np
import itertools

n = 6
iters = 1000
firstzero = 0
bothzero = 0
""" The next line iterates over arrays of length n+1 which contain only -1s and 1s """
for S in itertools.product([-1, 1], repeat=n+1):
    """For i from 0 to iters -1 """
    for i in xrange(iters):
        """ Choose a random array of length n.
            Prob 1/4 of being -1, prob 1/4 of being 1 and prob 1/2 of being 0. """
        F = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n)
        """The next loop just makes sure that F is not all zeros."""
        while np.all(F == 0):
            F = np.random.choice(np.array([-1, 0, 0, 1], dtype=np.int8), size=n)
        """np.convolve(F, S, 'valid') computes two inner products between
        F and the two successive windows of S of length n."""
        FS = np.convolve(F, S, 'valid')
        if FS[0] == 0:
            firstzero += 1
        if np.all(FS == 0):
            bothzero += 1

print("firstzero: %i" % firstzero)
print("bothzero: %i" % bothzero)

Ele está contando o número de vezes que a convolução de duas matrizes aleatórias, uma mais longa que a outra, com uma distribuição de probabilidade específica, tem um 0 na primeira posição ou um 0 em ambas as posições.

Eu fiz uma aposta com um amigo que diz que o Python é uma linguagem terrível para escrever código, que precisa ser rápido. Demora 9s no meu computador. Ele diz que pode ser feito 100 vezes mais rápido se escrito em uma "linguagem adequada".

O desafio é verificar se esse código pode ser feito 100 vezes mais rápido em qualquer idioma de sua escolha. Vou testar seu código e o mais rápido daqui a uma semana vencerá. Se alguém fica abaixo de 0,09s, ganha automaticamente e eu perco.

Status

  • Python . 30 vezes mais rápido por Alistair Buxon! Embora não seja a solução mais rápida, é de fato a minha favorita.
  • Oitava . 100 vezes mais rápido por @Thethos.
  • Ferrugem . 500 vezes mais rápido por @dbaupp.
  • C ++ . 570 vezes mais rápido por Guy Sirton.
  • C . 727 vezes mais rápido por @ace.
  • C ++ . Inacreditavelmente rápido por @Stefan.

As soluções mais rápidas agora são muito rápidas para um tempo sensato. Portanto, aumentei n para 10 e defina iters = 100000 para comparar os melhores. Sob essa medida, os mais rápidos são.

  • C . 7.5s por @ace.
  • C ++ . 1s por @Stefan.

Minha máquina Os horários serão executados na minha máquina. Esta é uma instalação padrão do ubuntu em um processador AMD FX-8350 de oito núcleos. Isso também significa que eu preciso poder executar seu código.

Acompanhamento publicado Como essa competição foi muito fácil de obter uma aceleração de x100, publiquei um acompanhamento para aqueles que desejam exercer seus conhecimentos em guru da velocidade. Veja o quão lento é realmente o Python (parte II)?

Comunidade
fonte

Respostas:

61

C ++ bit magic

0.84ms com RNG simples, 1.67ms com c ++ 11 std :: knuth

0.16ms com ligeira modificação algorítmica (veja editar abaixo)

A implementação do python é executada em 7,97 segundos no meu equipamento. Portanto, isso é 9488 a 4772 vezes mais rápido, dependendo do RNG que você escolher.

#include <iostream>
#include <bitset>
#include <random>
#include <chrono>
#include <stdint.h>
#include <cassert>
#include <tuple>

#if 0
// C++11 random
std::random_device rd;
std::knuth_b gen(rd());

uint32_t genRandom()
{
    return gen();
}
#else
// bad, fast, random.

uint32_t genRandom()
{
    static uint32_t seed = std::random_device()();
    auto oldSeed = seed;
    seed = seed*1664525UL + 1013904223UL; // numerical recipes, 32 bit
    return oldSeed;
}
#endif

#ifdef _MSC_VER
uint32_t popcnt( uint32_t x ){ return _mm_popcnt_u32(x); }
#else
uint32_t popcnt( uint32_t x ){ return __builtin_popcount(x); }
#endif



std::pair<unsigned, unsigned> convolve()
{
    const uint32_t n = 6;
    const uint32_t iters = 1000;
    unsigned firstZero = 0;
    unsigned bothZero = 0;

    uint32_t S = (1 << (n+1));
    // generate all possible N+1 bit strings
    // 1 = +1
    // 0 = -1
    while ( S-- )
    {
        uint32_t s1 = S % ( 1 << n );
        uint32_t s2 = (S >> 1) % ( 1 << n );
        uint32_t fmask = (1 << n) -1; fmask |= fmask << 16;
        static_assert( n < 16, "packing of F fails when n > 16.");


        for( unsigned i = 0; i < iters; i++ )
        {
            // generate random bit mess
            uint32_t F;
            do {
                F = genRandom() & fmask;
            } while ( 0 == ((F % (1 << n)) ^ (F >> 16 )) );

            // Assume F is an array with interleaved elements such that F[0] || F[16] is one element
            // here MSB(F) & ~LSB(F) returns 1 for all elements that are positive
            // and  ~MSB(F) & LSB(F) returns 1 for all elements that are negative
            // this results in the distribution ( -1, 0, 0, 1 )
            // to ease calculations we generate r = LSB(F) and l = MSB(F)

            uint32_t r = F % ( 1 << n );
            // modulo is required because the behaviour of the leftmost bit is implementation defined
            uint32_t l = ( F >> 16 ) % ( 1 << n );

            uint32_t posBits = l & ~r;
            uint32_t negBits = ~l & r;
            assert( (posBits & negBits) == 0 );

            // calculate which bits in the expression S * F evaluate to +1
            unsigned firstPosBits = ((s1 & posBits) | (~s1 & negBits));
            // idem for -1
            unsigned firstNegBits = ((~s1 & posBits) | (s1 & negBits));

            if ( popcnt( firstPosBits ) == popcnt( firstNegBits ) )
            {
                firstZero++;

                unsigned secondPosBits = ((s2 & posBits) | (~s2 & negBits));
                unsigned secondNegBits = ((~s2 & posBits) | (s2 & negBits));

                if ( popcnt( secondPosBits ) == popcnt( secondNegBits ) )
                {
                    bothZero++;
                }
            }
        }
    }

    return std::make_pair(firstZero, bothZero);
}

int main()
{
    typedef std::chrono::high_resolution_clock clock;
    int rounds = 1000;
    std::vector< std::pair<unsigned, unsigned> > out(rounds);

    // do 100 rounds to get the cpu up to speed..
    for( int i = 0; i < 10000; i++ )
    {
        convolve();
    }


    auto start = clock::now();

    for( int i = 0; i < rounds; i++ )
    {
        out[i] = convolve();
    }

    auto end = clock::now();
    double seconds = std::chrono::duration_cast< std::chrono::microseconds >( end - start ).count() / 1000000.0;

#if 0
    for( auto pair : out )
        std::cout << pair.first << ", " << pair.second << std::endl;
#endif

    std::cout << seconds/rounds*1000 << " msec/round" << std::endl;

    return 0;
}

Compile em 64 bits para registros extras. Ao usar o gerador aleatório simples, os loops em convolve () são executados sem nenhum acesso à memória, todas as variáveis ​​são armazenadas nos registradores.

Como funciona: em vez de armazenar Se Fcomo matrizes na memória, ele é armazenado como bits em um uint32_t.
Pois S, os nbits menos significativos são usados ​​onde um bit definido indica um +1 e um bit não definido indica -1.
Frequer pelo menos 2 bits para criar uma distribuição de [-1, 0, 0, 1]. Isso é feito gerando bits aleatórios e examinando os 16 bits menos significativos (chamados r) e os 16 bits mais significativos (chamados l). Se l & ~rassumirmos que F é +1, se ~l & rassumirmos que Fé -1. Caso contrário, Fé 0. Isso gera a distribuição que estamos procurando.

Agora temos S, posBitscom um bit definido em todos os locais onde F == 1 e negBitscom um bit definido em todos os locais onde F == -1.

Podemos provar que F * S(onde * denota multiplicação) é avaliado como +1 sob a condição (S & posBits) | (~S & negBits). Também podemos gerar lógica semelhante para todos os casos em que é F * Savaliado como -1. E, finalmente, sabemos quesum(F * S) avalia como 0 se e somente se houver uma quantidade igual de -1 e + 1 no resultado. É muito fácil calcular isso simplesmente comparando o número de +1 e -1 bits.

Esta implementação usa 32 bits de ints, e o valor máximo n aceito é 16. É possível escalar a implementação para 31 bits modificando o código de geração aleatória e para 63 bits usando uint64_t em vez de uint32_t.

editar

A seguinte função convolve:

std::pair<unsigned, unsigned> convolve()
{
    const uint32_t n = 6;
    const uint32_t iters = 1000;
    unsigned firstZero = 0;
    unsigned bothZero = 0;
    uint32_t fmask = (1 << n) -1; fmask |= fmask << 16;
    static_assert( n < 16, "packing of F fails when n > 16.");


    for( unsigned i = 0; i < iters; i++ )
    {
        // generate random bit mess
        uint32_t F;
        do {
            F = genRandom() & fmask;
        } while ( 0 == ((F % (1 << n)) ^ (F >> 16 )) );

        // Assume F is an array with interleaved elements such that F[0] || F[16] is one element
        // here MSB(F) & ~LSB(F) returns 1 for all elements that are positive
        // and  ~MSB(F) & LSB(F) returns 1 for all elements that are negative
        // this results in the distribution ( -1, 0, 0, 1 )
        // to ease calculations we generate r = LSB(F) and l = MSB(F)

        uint32_t r = F % ( 1 << n );
        // modulo is required because the behaviour of the leftmost bit is implementation defined
        uint32_t l = ( F >> 16 ) % ( 1 << n );

        uint32_t posBits = l & ~r;
        uint32_t negBits = ~l & r;
        assert( (posBits & negBits) == 0 );

        uint32_t mask = posBits | negBits;
        uint32_t totalBits = popcnt( mask );
        // if the amount of -1 and +1's is uneven, sum(S*F) cannot possibly evaluate to 0
        if ( totalBits & 1 )
            continue;

        uint32_t adjF = posBits & ~negBits;
        uint32_t desiredBits = totalBits / 2;

        uint32_t S = (1 << (n+1));
        // generate all possible N+1 bit strings
        // 1 = +1
        // 0 = -1
        while ( S-- )
        {
            // calculate which bits in the expression S * F evaluate to +1
            auto firstBits = (S & mask) ^ adjF;
            auto secondBits = (S & ( mask << 1 ) ) ^ ( adjF << 1 );

            bool a = desiredBits == popcnt( firstBits );
            bool b = desiredBits == popcnt( secondBits );
            firstZero += a;
            bothZero += a & b;
        }
    }

    return std::make_pair(firstZero, bothZero);
}

reduz o tempo de execução para 0,160-0,161ms. O desenrolar manual do loop (não na foto acima) faz com que 0,150. O caso menos trivial de n = 10, iter = 100000 é executado em 250ms. Tenho certeza de que posso obtê-lo abaixo de 50ms usando núcleos adicionais, mas isso é muito fácil.

Isso é feito liberando o ramo do loop interno e trocando o loop F e S.
Se bothZeronão for necessário, eu posso reduzir o tempo de execução para 0,02 ms, fazendo um loop esparso sobre todas as matrizes S possíveis.

Stefan
fonte
3
Você poderia fornecer uma versão amigável para o gcc e também qual seria sua linha de comando? Não tenho certeza se posso testá-lo atualmente.
Não sei nada sobre isso, mas o Google me diz que __builtin_popcount pode ser um substituto para _mm_popcnt_u32 ().
3
Código atualizado, usa a opção #ifdef para selecionar o comando popcnt correto. Ele compila -std=c++0x -mpopcnt -O2e leva 1.01ms para executar no modo de 32 bits (não tenho uma versão do GCC de 64 bits em mãos).
Stefan
Você poderia imprimir a saída? Eu não tenho certeza se ele está realmente fazendo atualmente qualquer coisa :)
7
Você é claramente um mago. + 1
BurntPizza
76

Python2.7 + Numpy 1.8.1: 10.242 s

Fortran 90+: 0,029 s 0,003 s 0,022 s 0,010 s

Droga, você perdeu sua aposta! Também não é uma gota de paralelismo, apenas o Fortran 90+ direto.

Edição Eu peguei o algoritmo de Guy Sirton para permutar o array S(boa descoberta: D). Aparentemente, eu também tinha os -g -tracebacksinalizadores do compilador ativos, que estavam diminuindo esse código para cerca de 0,017s. Atualmente, estou compilando isso como

ifort -fast -o convolve convolve_random_arrays.f90

Para quem não tem ifort, você pode usar

gfortran -O3 -ffast-math -o convolve convolve_random_arrays.f90

EDIT 2 : A diminuição no tempo de execução é porque eu estava fazendo algo errado anteriormente e obtive uma resposta incorreta. Fazer da maneira certa é aparentemente mais lento. Ainda não consigo acreditar que o C ++ seja mais rápido que o meu, por isso provavelmente vou gastar algum tempo esta semana tentando ajustar essa porcaria para acelerar.

EDIÇÃO 3 : Simplesmente alterando a seção RNG usando uma baseada no RNG da BSD (como sugerido por Sampo Smolander) e eliminando a divisão constante por m1, reduzi o tempo de execução para o mesmo que a resposta C ++ de Guy Sirton . O uso de matrizes estáticas (como sugerido pela Sharpie) reduz o tempo de execução para o tempo de execução do C ++! Yay Fortran! : D

EDIT 4 Aparentemente, isso não é compilado (com gfortran) e executado corretamente (valores incorretos) porque os números inteiros estão ultrapassando seus limites. Fiz correções para garantir que funcione, mas isso requer que se tenha ifort 11+ ou gfortran 4.7+ (ou outro compilador que permita iso_fortran_enve o int64tipo F2008 ).

Aqui está o código:

program convolve_random_arrays
   use iso_fortran_env
   implicit none
   integer(int64), parameter :: a1 = 1103515245
   integer(int64), parameter :: c1 = 12345
   integer(int64), parameter :: m1 = 2147483648
   real, parameter ::    mi = 4.656612873e-10 ! 1/m1
   integer, parameter :: n = 6
   integer :: p, pmax, iters, i, nil(0:1), seed
   !integer, allocatable ::  F(:), S(:), FS(:)
   integer :: F(n), S(n+1), FS(2)

   !n = 6
   !allocate(F(n), S(n+1), FS(2))
   iters = 1000
   nil = 0

   !call init_random_seed()

   S = -1
   pmax = 2**(n+1)
   do p=1,pmax
      do i=1,iters
         F = rand_int_array(n)
         if(all(F==0)) then
            do while(all(F==0))
               F = rand_int_array(n)
            enddo
         endif

         FS = convolve(F,S)

         if(FS(1) == 0) then
            nil(0) = nil(0) + 1
            if(FS(2) == 0) nil(1) = nil(1) + 1
         endif

      enddo
      call permute(S)
   enddo

   print *,"first zero:",nil(0)
   print *," both zero:",nil(1)

 contains
   pure function convolve(x, h) result(y)
!x is the signal array
!h is the noise/impulse array
      integer, dimension(:), intent(in) :: x, h
      integer, dimension(abs(size(x)-size(h))+1) :: y
      integer:: i, j, r
      y(1) = dot_product(x,h(1:n-1))
      y(2) = dot_product(x,h(2:n  ))
   end function convolve

   pure subroutine permute(x)
      integer, intent(inout) :: x(:)
      integer :: i

      do i=1,size(x)
         if(x(i)==-1) then
            x(i) = 1
            return
         endif
         x(i) = -1
      enddo
   end subroutine permute

   function rand_int_array(i) result(x)
     integer, intent(in) :: i
     integer :: x(i), j
     real :: y
     do j=1,i
        y = bsd_rng()
        if(y <= 0.25) then
           x(j) = -1
        else if (y >= 0.75) then
           x(j) = +1
        else
           x(j) = 0
        endif
     enddo
   end function rand_int_array

   function bsd_rng() result(x)
      real :: x
      integer(int64) :: b=3141592653
      b = mod(a1*b + c1, m1)
      x = real(b)*mi
   end function bsd_rng
end program convolve_random_arrays

Suponho que a pergunta agora é: você irá parar de usar o Python lento como melaço e usar o Fortran rápido como elétrons que pode se mover;).

Kyle Kanos
fonte
1
A declaração de caso não seria mais rápida do que uma função de gerador? A menos que você esteja esperando algum tipo de aceleração de previsão de ramificação / linha de cache / etc?
precisa saber é o seguinte
17
A velocidade deve ser comparada na mesma máquina. Qual tempo de execução você obteve para o código do OP?
Nbubis
3
A resposta C ++ implementa seu próprio gerador de números aleatórios muito leve. Sua resposta usou o padrão que acompanha o compilador, o que poderia ser mais lento?
Sampo Smolander
3
Além disso, o exemplo C ++ parece estar usando matrizes alocadas estaticamente. Tente usar matrizes de comprimento fixo definidas em tempo de compilação e veja se ela se depila em algum momento.
Sharpie
1
@KyleKanos @Lembik o problema é que a atribuição de número inteiro no fortran não está usando implicitamente a especificação int64, portanto, os números são int32 antes de qualquer conversão. O código deve ser: integer(int64) :: b = 3141592653_int64para todos os int64. Isso faz parte do padrão fortran e é esperado pelo programador em uma linguagem de programação declarada por tipo. (Note que as configurações padrão de curso pode substituir esse)
zeroth
69

Python 2.7 - 0.882s 0.283s

(Original do OP: 6.404s)

Edit: Otimização de Steven Rumbalski pré-computando valores de F. Com essa otimização, o cpython supera os 0.365s do pypy.

import itertools
import operator
import random

n=6
iters = 1000
firstzero = 0
bothzero = 0

choicesF = filter(any, itertools.product([-1, 0, 0, 1], repeat=n))

for S in itertools.product([-1,1], repeat = n+1):
    for i in xrange(iters):
        F = random.choice(choicesF)
        if not sum(map(operator.mul, F, S[:-1])):
            firstzero += 1
            if not sum(map(operator.mul, F, S[1:])):
                bothzero += 1

print "firstzero", firstzero
print "bothzero", bothzero

O código original do OP usa matrizes tão pequenas que não há benefício em usar o Numpy, como demonstra a implementação pura do python. Mas veja também essa implementação numpy que é três vezes mais rápida do que meu código.

Também otimizo pulando o resto da convolução se o primeiro resultado não for zero.

Alistair Buxton
fonte
11
Com o pypy, isso é executado em cerca de 0,5 segundos.
Alistair Buxton
2
Você obtém uma aceleração muito mais convincente se definir n = 10. Recebo 19s versus 4,6s para cpython versus pypy.
3
Outra otimização seria precomputar as possibilidades Fporque existem apenas 4032 delas. Defina choicesF = filter(any, itertools.product([-1, 0, 0, 1], repeat=n))fora dos loops. Em seguida, no loop interno, defina F = random.choice(choicesF). Eu recebo uma aceleração de 3x com essa abordagem.
precisa saber é o seguinte
3
Que tal compilar isso no Cython? Então, adicionando alguns tipos estáticos de tato?
Thane Brimhall
2
Coloque tudo em uma função e chame no final. Isso localiza os nomes, o que também faz com que a otimização proposta pelo @riffraff funcione. Além disso, mova a criação de range(iters)fora do loop. Ao todo, recebo uma aceleração de cerca de 7% em relação à sua ótima resposta.
WolframH
44

Ferrugem: 0.011s

Python original: 8.3

Uma tradução direta do Python original.

extern crate rand;

use rand::Rng;

static N: uint = 6;
static ITERS: uint = 1000;

fn convolve<T: Num>(into: &mut [T], a: &[T], b: &[T]) {
    // we want `a` to be the longest array
    if a.len() < b.len() {
        convolve(into, b, a);
        return
    }

    assert_eq!(into.len(), a.len() - b.len() + 1);

    for (n,place) in into.mut_iter().enumerate() {
        for (x, y) in a.slice_from(n).iter().zip(b.iter()) {
            *place = *place + *x * *y
        }
    }
}

fn main() {
    let mut first_zero = 0;
    let mut both_zero = 0;
    let mut rng = rand::XorShiftRng::new().unwrap();

    for s in PlusMinus::new() {
        for _ in range(0, ITERS) {
            let mut f = [0, .. N];
            while f.iter().all(|x| *x == 0) {
                for p in f.mut_iter() {
                    match rng.gen::<u32>() % 4 {
                        0 => *p = -1,
                        1 | 2 => *p = 0,
                        _ => *p = 1
                    }
                }
            }

            let mut fs = [0, .. 2];
            convolve(fs, s, f);

            if fs[0] == 0 { first_zero += 1 }
            if fs.iter().all(|&x| x == 0) { both_zero += 1 }
        }
    }

    println!("{}\n{}", first_zero, both_zero);
}



/// An iterator over [+-]1 arrays of the appropriate length
struct PlusMinus {
    done: bool,
    current: [i32, .. N + 1]
}
impl PlusMinus {
    fn new() -> PlusMinus {
        PlusMinus { done: false, current: [-1, .. N + 1] }
    }
}

impl Iterator<[i32, .. N + 1]> for PlusMinus {
    fn next(&mut self) -> Option<[i32, .. N+1]> {
        if self.done {
            return None
        }

        let ret = self.current;

        // a binary "adder", that just adds one to a bit vector (where
        // -1 is the zero, and 1 is the one).
        for (i, place) in self.current.mut_iter().enumerate() {
            *place = -*place;
            if *place == 1 {
                break
            } else if i == N {
                // we've wrapped, so we want to stop after this one
                self.done = true
            }
        }

        Some(ret)
    }
}
  • Compilado com --opt-level=3
  • Meu compilador de ferrugem é recente todas as noites : ( rustc 0.11-pre-nightly (eea4909 2014-04-24 23:41:15 -0700)para ser mais preciso)
huon
fonte
Eu consegui compilar usando a versão noturna da ferrugem. No entanto, acho que o código está errado. A saída deve ser algo próximo a firstzero 27215 bothzero 12086. Em vez disso, dá 27367 6481
@Lembik, whoops, misturou meus ae bs na convolução; fixo (não altera notavelmente o tempo de execução).
huon
4
É uma demonstração muito boa da velocidade da ferrugem.
39

C ++ (VS 2012) - 0.026s 0.015s

Python 2.7.6 / Numpy 1.8.1 - 12s

Aceleração ~ x800.

A diferença seria muito menor se as matrizes envolvidas fossem muito grandes ...

#include <vector>
#include <iostream>
#include <ctime>

using namespace std;

static unsigned int seed = 35;

int my_random()
{
   seed = seed*1664525UL + 1013904223UL; // numerical recipes, 32 bit

   switch((seed>>30) & 3)
   {
   case 0: return 0;
   case 1: return -1;
   case 2: return 1;
   case 3: return 0;
   }
   return 0;
}

bool allzero(const vector<int>& T)
{
   for(auto x : T)
   {
      if(x!=0)
      {
         return false;
      }
   }
   return true;
}

void convolve(vector<int>& out, const vector<int>& v1, const vector<int>& v2)
{
   for(size_t i = 0; i<out.size(); ++i)
   {
      int result = 0;
      for(size_t j = 0; j<v2.size(); ++j)
      {
         result += v1[i+j]*v2[j];
      }
      out[i] = result;
   }
}

void advance(vector<int>& v)
{
   for(auto &x : v)
   {
      if(x==-1)
      {
         x = 1;
         return;
      }
      x = -1;
   }
}

void convolve_random_arrays(void)
{
   const size_t n = 6;
   const int two_to_n_plus_one = 128;
   const int iters = 1000;
   int bothzero = 0;
   int firstzero = 0;

   vector<int> S(n+1);
   vector<int> F(n);
   vector<int> FS(2);

   time_t current_time;
   time(&current_time);
   seed = current_time;

   for(auto &x : S)
   {
      x = -1;
   }
   for(int i=0; i<two_to_n_plus_one; ++i)
   {
      for(int j=0; j<iters; ++j)
      {
         do
         {
            for(auto &x : F)
            {
               x = my_random();
            }
         } while(allzero(F));
         convolve(FS, S, F);
         if(FS[0] == 0)
         {
            firstzero++;
            if(FS[1] == 0)
            {
               bothzero++;
            }
         }
      }
      advance(S);
   }
   cout << firstzero << endl; // This output can slow things down
   cout << bothzero << endl; // comment out for timing the algorithm
}

Algumas notas:

  • A função aleatória está sendo chamada no loop, então procurei um gerador congruencial linear muito leve (mas olhei generosamente para os MSBs).
  • Este é realmente apenas o ponto de partida para uma solução otimizada.
  • Não demorou tanto tempo para escrever ...
  • Eu percorro todos os valores de S considerando S[0]o dígito "menos significativo".

Adicione esta função principal para um exemplo independente:

int main(int argc, char** argv)
{
  for(int i=0; i<1000; ++i) // run 1000 times for stop-watch
  {
      convolve_random_arrays();
  }
}
Guy Sirton
fonte
1
De fato. O tamanho minúsculo das matrizes no código do OP significa que usar numpy é na verdade uma ordem de magnitude mais lenta que o python direto.
Alistair Buxton
2
Agora x800 é o que eu estou falando!
Muito agradável! Eu já aumentou a velocidade-up no meu código por causa de sua advancefunção, assim meu código é agora mais rápido do que o seu: P (mas muito boa competição!)
Kyle Kanos
1
@embemb sim, como diz Mat. Você precisa de C ++ 11 supprt e uma função principal. Deixe-me saber se você precisar de mais ajuda para obter este a correr ...
Guy Sirton
2
Eu só testei isso e poderia raspar de outros 20% usando matrizes simples em vez de std :: vector ..
PlasmaHH
21

C

Leva 0,015s na minha máquina, com o código original do OP demorando ~ 7,7s. Tentei otimizar gerando a matriz aleatória e convolvendo no mesmo loop, mas isso não parece fazer muita diferença.

A primeira matriz é gerada usando um número inteiro, escreva-o em binário e altere todos os 1 para -1 e todos os 0 para 1. O restante deve ser bem direto.

Edit: em vez de ter ncomo uma int, agora temos ncomo uma constante definida por macro, para que possamos usar em int arr[n];vez de malloc.

Edit2: Em vez de rand()função interna, isso agora implementa um PRNG xorshift. Além disso, muitas instruções condicionais são removidas ao gerar a matriz aleatória.

Instruções de compilação:

gcc -O3 -march=native -fwhole-program -fstrict-aliasing -ftree-vectorize -Wall ./test.c -o ./test

Código:

#include <stdio.h>
#include <time.h>

#define n (6)
#define iters (1000)
unsigned int x,y=34353,z=57768,w=1564; //PRNG seeds

/* xorshift PRNG
 * Taken from https://en.wikipedia.org/wiki/Xorshift#Example_implementation
 * Used under CC-By-SA */
int myRand() {
    unsigned int t;
    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = w ^ (w >> 19) ^ t ^ (t >> 8);
}

int main() {
    int firstzero=0, bothzero=0;
    int arr[n+1];
    unsigned int i, j;
    x=(int)time(NULL);

    for(i=0; i< 1<<(n+1) ; i++) {
        unsigned int tmp=i;
        for(j=0; j<n+1; j++) {
            arr[j]=(tmp&1)*(-2)+1;
            tmp>>=1;
        }
        for(j=0; j<iters; j++) {
            int randArr[n];
            unsigned int k, flag=0;
            int first=0, second=0;
            do {
                for(k=0; k<n; k++) {
                    randArr[k]=(1-(myRand()&3))%2;
                    flag+=(randArr[k]&1);
                    first+=arr[k]*randArr[k];
                    second+=arr[k+1]*randArr[k];
                }
            } while(!flag);
            firstzero+=(!first);
            bothzero+=(!first&&!second);
        }
    }
    printf("firstzero %d\nbothzero %d\n", firstzero, bothzero);
    return 0;
}
ace_HongKongIndependence
fonte
1
Eu testei isso. É muito rápido (tente n = 10) e fornece saída com aparência correta. Obrigado.
Essa implementação não segue o original porque, se o vetor aleatório tiver todos os zeros, apenas o último elemento seria gerado novamente. No original, o vetor inteiro seria. Você precisa anexar esse loop do{}while(!flag)ou algo nesse sentido. Não espero que mude muito o tempo de execução (pode torná-lo mais rápido).
precisa saber é o seguinte
@Guy Sirton Note que antes da continue;declaração que eu atribuído -1a k, então kfará um loop de 0 novamente.
ace_HongKongIndependence
1
@ace ah! você está certo. Eu estava digitalizando muito rápido e parecia que era -=melhor que =-:-) Um loop while seria mais legível.
precisa saber é o seguinte
17

J

Eu não espero derrotar nenhum idioma compilado, e algo me diz que seria necessário uma máquina milagrosa para obter menos de 0,09 s com isso, mas eu gostaria de enviar esse J de qualquer maneira, porque é muito liso.

NB. constants
num =: 6
iters =: 1000

NB. convolve
NB. take the multiplication table                */
NB. then sum along the NE-SW diagonals           +//.
NB. and keep the longest ones                    #~ [: (= >./) #/.
NB. operate on rows of higher dimensional lists  " 1
conv =: (+//. #~ [: (= >./) #/.) @: (*/) " 1

NB. main program
S  =: > , { (num+1) # < _1 1                NB. all {-1,1}^(num+1)
F  =: (3&= - 0&=) (iters , num) ?@$ 4       NB. iters random arrays of length num
FS =: ,/ S conv/ F                          NB. make a convolution table
FB =: +/ ({. , *./)"1 ] 0 = FS              NB. first and both zero
('first zero ',:'both zero ') ,. ":"0 FB    NB. output results

Isso leva cerca de 0,5 s em um laptop da década anterior, apenas 20x mais rápido que o Python na resposta. A maior parte do tempo é gasta convporque a escrevemos preguiçosamente (calculamos toda a convolução) e com total generalidade.

Como sabemos Se Fpodemos acelerar, otimizações específicas para este programa. O melhor que pude apresentar é - conv =: ((num, num+1) { +//.)@:(*/)"1selecione especificamente os dois números que correspondem das somas diagonais aos elementos mais longos da convolução - que reduzem pela metade o tempo.

algoritmshark
fonte
6
Sempre vale a pena enviar J, cara :) #
Vitaly Dyatlov
17

Perl - 9.3X mais rápido ... 830% de melhoria

No meu netbook antigo, o código do OP demora 53 segundos para ser executado; A versão de Alistair Buxton leva cerca de 6,5 segundos e a versão Perl a seguir leva cerca de 5,7 segundos.

use v5.10;
use strict;
use warnings;

use Algorithm::Combinatorics qw( variations_with_repetition );
use List::Util qw( any sum );
use List::MoreUtils qw( pairwise );

my $n         = 6;
my $iters     = 1000;
my $firstzero = 0;
my $bothzero  = 0;

my $variations = variations_with_repetition([-1, 1], $n+1);
while (my $S = $variations->next)
{
  for my $i (1 .. $iters)
  {
    my @F;
    until (@F and any { $_ } @F)
    {
      @F = map +((-1,0,0,1)[rand 4]), 1..$n;
    }

    # The pairwise function doesn't accept array slices,
    # so need to copy into a temp array @S0
    my @S0 = @$S[0..$n-1];

    unless (sum pairwise { $a * $b } @F, @S0)
    {
      $firstzero++;
      my @S1 = @$S[1..$n];  # copy again :-(
      $bothzero++ unless sum pairwise { $a * $b } @F, @S1;
    }
  }
}

say "firstzero ", $firstzero;
say "bothzero ", $bothzero;
tobyink
fonte
12

Python 2.7 - numpy 1.8.1 com ligações mkl - 0.086s

(Original do OP: 6.404s) (python puro de Buxton: 0.270s)

import numpy as np
import itertools

n=6
iters = 1000

#Pack all of the Ses into a single array
S = np.array( list(itertools.product([-1,1], repeat=n+1)) )

# Create a whole array of test arrays, oversample a bit to ensure we 
# have at least (iters) of them
F = np.random.rand(int(iters*1.1),n)
F = ( F < 0.25 )*-1 + ( F > 0.75 )*1
goodrows = (np.abs(F).sum(1)!=0)
assert goodrows.sum() > iters, "Got very unlucky"
# get 1000 cases that aren't all zero
F = F[goodrows][:iters]

# Do the convolution explicitly for the two 
# slots, but on all of the Ses and Fes at the 
# same time
firstzeros = (F[:,None,:]*S[None,:,:-1]).sum(-1)==0
secondzeros = (F[:,None,:]*S[None,:,1:]).sum(-1)==0

firstzero_count = firstzeros.sum()
bothzero_count = (firstzeros * secondzeros).sum()
print "firstzero", firstzero_count
print "bothzero", bothzero_count

Como Buxton aponta, o código original do OP usa matrizes tão pequenas que não há benefício em usar o Numpy. Essa implementação aproveita o numpy, executando todos os casos F e S de uma só vez, de maneira orientada por matriz. Isso combinado com as ligações mkl para python leva a uma implementação muito rápida.

Observe também que apenas carregar as bibliotecas e iniciar o intérprete leva 0,076s; portanto, o cálculo real leva ~ 0,01 segundos, semelhante à solução C ++.

alemi
fonte
O que são ligações mkl e como as consigo no ubuntu?
Execução python -c "import numpy; numpy.show_config()"irá mostrar-lhe se a sua versão do numpy é compilado contra blas / atlas / mkl, etc. ATLAS é um pacote de matemática acelerada livre que numpy pode ser ligado contra , Intel MKL você geralmente tem que pagar (a menos que você é um acadêmico) e pode ser vinculado a numpy / scipy .
alemi
Para facilitar, use a distribuição anaconda python e use o pacote accelerate . Ou use a distribuição entusiasmada .
alemi
Se você estiver no Windows, faça o download do numpy aqui . Instaladores numpy pré-compilados vinculados ao MKL.
Nome falso
9

MATLAB 0.024s

Computador 1

  • Código original: ~ 3.3 s
  • Código de Alistar Buxton: ~ 0.51 s
  • Novo código de Alistar Buxton: ~ 0,25 s
  • Código Matlab: ~ 0.024 s (o Matlab já está em execução)

Computador 2

  • Código original: ~ 6.66 s
  • Código de Alistar Buxton: ~ 0.64 s
  • Novo código de Alistar Buxton:?
  • Matlab: ~ 0,07 s (o Matlab já está em execução)
  • Oitava: ~ 0.07 s

Decidi experimentar o Matlab tão lento. Se você sabe como, pode se livrar da maioria dos loops (no Matlab), o que o torna muito rápido. No entanto, os requisitos de memória são mais altos do que as soluções em loop, mas isso não será um problema se você não tiver matrizes muito grandes ...

function call_convolve_random_arrays
tic
convolve_random_arrays
toc
end

function convolve_random_arrays

n = 6;
iters = 1000;
firstzero = 0;
bothzero = 0;

rnd = [-1, 0, 0, 1];

S = -1 *ones(1, n + 1);

IDX1 = 1:n;
IDX2 = IDX1 + 1;

for i = 1:2^(n + 1)
    F = rnd(randi(4, [iters, n]));
    sel = ~any(F,2);
    while any(sel)
        F(sel, :) = rnd(randi(4, [sum(sel), n]));
        sel = ~any(F,2);
    end

    sum1 = F * S(IDX1)';
    sel = sum1 == 0;
    firstzero = firstzero + sum(sel);

    sum2 = F(sel, :) * S(IDX2)';
    sel = sum2 == 0;
    bothzero = bothzero + sum(sel);

    S = permute(S); 
end

fprintf('firstzero %i \nbothzero %i \n', firstzero, bothzero);

end

function x = permute(x)

for i=1:length(x)
    if(x(i)==-1)
        x(i) = 1;
            return
    end
        x(i) = -1;
end

end

Aqui está o que eu faço:

  • use a função Kyle Kanos para permutar através de S
  • calcular todos os números aleatórios n * iters de uma só vez
  • mapa 1 a 4 a [-1 0 0 1]
  • use a multiplicação de matrizes (soma dos elementos (F * S (1: 5)) é igual à multiplicação de matrizes de F * S (1: 5) '
  • para bothzero: calcule apenas membros que cumpram a primeira condição

Presumo que você não tenha o Matlab, o que é muito ruim, pois eu realmente gostaria de ver como ele se compara ...

(A função pode ser mais lenta na primeira vez que você a executa.)

Mathause
fonte
Bem, eu tenho oitava se você pode fazê-lo funcionar para isso ...?
Eu posso tentar - eu nunca trabalhei com oitava, no entanto.
mathause
Ok, eu posso executá-lo como está na oitava se eu colocar o código em um arquivo chamado call_convolve_random_arrays.m e depois chamá-lo da oitava.
mathause
Precisa de mais código para realmente fazer alguma coisa? Quando eu faço "octave call_convolve_random_arrays.m", ele não gera nada. Veja bpaste.net/show/JPtLOCeI3aP3wc3F3aGf
desculpe, tente abrir a oitava e execute-a então. Ele deve exibir firstzero, bothzero e tempo de execução.
mathause
7

Julia: 0,30 s

Op's Python: 21.36 s (Core2 duo)

71x aceleração

function countconv()                                                                                                                                                           
    n = 6                                                                                                                                                                      
    iters = 1000                                                                                                                                                               
    firstzero = 0                                                                                                                                                              
    bothzero = 0                                                                                                                                                               
    cprod= Iterators.product(fill([-1,1], n+1)...)                                                                                                                             
    F=Array(Float64,n);                                                                                                                                                        
    P=[-1. 0. 0. 1.]                                                                                                                                                                                                                                                                                                             

    for S in cprod                                                                                                                                                             
        Sm=[S...]                                                                                                                                                              
        for i = 1:iters                                                                                                                                                        
            F=P[rand(1:4,n)]                                                                                                                                                  
            while all(F==0)                                                                                                                                                   
                F=P[rand(1:4,n)]                                                                                                                                              
            end                                                                                                                                                               
            if  dot(reverse!(F),Sm[1:end-1]) == 0                                                                                                                           
                firstzero += 1                                                                                                                                                 
                if dot(F,Sm[2:end]) == 0                                                                                                                              
                    bothzero += 1                                                                                                                                              
                end                                                                                                                                                            
            end                                                                                                                                                                
        end                                                                                                                                                                    
    end
    return firstzero,bothzero
end

Fiz algumas modificações na resposta de Arman na Julia: Antes de tudo, envolvi-a em uma função, pois as variáveis ​​globais dificultam a inferência de tipo de Julia e o JIT: uma variável global pode alterar seu tipo a qualquer momento e deve ser verificada a cada operação . Então, me livrei das funções anônimas e da compreensão de arrays. Eles não são realmente necessários e ainda são muito lentos. Julia é mais rápida com abstrações de nível inferior no momento.

Existem muitas outras maneiras de torná-lo mais rápido, mas isso faz um trabalho decente.

user20768
fonte
Você está medindo o tempo no REPL ou executando o arquivo inteiro na linha de comando?
Aditya
ambos do REPL.
user20768
6

Ok, estou postando isso apenas porque sinto que o Java precisa ser representado aqui. Sou péssimo com outras línguas e confesso que não entendi exatamente o problema, por isso precisarei de ajuda para corrigir esse código. Roubei a maior parte do exemplo C do ás de código e depois peguei emprestados alguns trechos de outros. Espero que não seja um falso passo ...

Uma coisa que eu gostaria de destacar é que os idiomas que otimizam em tempo de execução precisam ser executados várias / muitas vezes para atingir a velocidade máxima. Eu acho que é justificável usar a velocidade totalmente otimizada (ou pelo menos a velocidade média) porque a maioria das coisas com as quais você está preocupado em correr rápido será executada várias vezes.

O código ainda precisa ser corrigido, mas eu o executei de qualquer maneira para ver que horas chegaria.

Aqui estão os resultados de uma CPU Intel (R) Xeon (E) E3-1270 V2 a 3.50GHz no Ubuntu executando-a 1000 vezes:

servidor: / tmp # time java8 -cp. Testador

firstzero 40000

bothzero 20000

primeiro tempo de execução: 41 ms último tempo de execução: 4 ms

usuário 0m5.014s real 0m4.664s sys 0m0.268s

Aqui está o meu código de baixa qualidade:

public class Tester 
{
    public static void main( String[] args )
    {
        long firstRunTime = 0;
        long lastRunTime = 0;
        String testResults = null;
        for( int i=0 ; i<1000 ; i++ )
        {
            long timer = System.currentTimeMillis();
            testResults = new Tester().runtest();
            lastRunTime = System.currentTimeMillis() - timer;
            if( i ==0 )
            {
                firstRunTime = lastRunTime;
            }
        }
        System.err.println( testResults );
        System.err.println( "first run time: " + firstRunTime + " ms" );
        System.err.println( "last run time: " + lastRunTime + " ms" );
    }

    private int x,y=34353,z=57768,w=1564; 

    public String runtest()
    {
        int n = 6;
        int iters = 1000;
        //#define iters (1000)
        //PRNG seeds

        /* xorshift PRNG
         * Taken from https://en.wikipedia.org/wiki/Xorshift#Example_implementation
         * Used under CC-By-SA */

            int firstzero=0, bothzero=0;
            int[] arr = new int[n+1];
            int i=0, j=0;
            x=(int)(System.currentTimeMillis()/1000l);

            for(i=0; i< 1<<(n+1) ; i++) {
                int tmp=i;
                for(j=0; j<n+1; j++) {
                    arr[j]=(tmp&1)*(-2)+1;
                    tmp>>=1;
                }
                for(j=0; j<iters; j++) {
                    int[] randArr = new int[n];
                    int k=0;
                    long flag = 0;
                    int first=0, second=0;
                    do {
                        for(k=0; k<n; k++) {
                            randArr[k]=(1-(myRand()&3))%2;
                            flag+=(randArr[k]&1);
                            first+=arr[k]*randArr[k];
                            second+=arr[k+1]*randArr[k];
                        }
                    } while(allzero(randArr));
                    if( first == 0 )
                    {
                        firstzero+=1;
                        if( second == 0 )
                        {
                            bothzero++;
                        }
                    }
                }
            }
         return ( "firstzero " + firstzero + "\nbothzero " + bothzero + "\n" );
    }

    private boolean allzero(int[] arr)
    {
       for(int x : arr)
       {
          if(x!=0)
          {
             return false;
          }
       }
       return true;
    }

    public int myRand() 
    {
        long t;
        t = x ^ (x << 11);
        x = y; y = z; z = w;
        return (int)( w ^ (w >> 19) ^ t ^ (t >> 8));
    }
}

E tentei executar o código python depois de atualizar o python e instalar o python-numpy, mas recebo o seguinte:

server:/tmp# python tester.py
Traceback (most recent call last):
  File "peepee.py", line 15, in <module>
    F = np.random.choice(np.array([-1,0,0,1], dtype=np.int8), size = n)
AttributeError: 'module' object has no attribute 'choice'
Chris Seline
fonte
Comentários: Nunca use currentTimeMillispara benchmarking (use a versão nano no System) e as execuções de 1k podem não ser suficientes para envolver o JIT (1,5k para o cliente e 10k para o servidor seriam os padrões, embora você chame o myRand com freqüência suficiente JITed que deve fazer com que algumas funções do callstack sejam compiladas, o que pode funcionar aqui). Última, mas não menos importante, o PNRG fraco está trapaceando, mas o mesmo acontece com a solução C ++ e outras, então acho que isso não é injusto demais.
Voo
No Windows, você precisa evitar o currentTimeMillis, mas para o linux, para todas as medições de granularidade, exceto muito finas, você não precisa de tempo nano, e a chamada para obter tempo nano é muito mais cara que milis. Então, eu discordo muito que você nunca deve usá-lo.
Chris Seline 5/05
Então, você está escrevendo código Java para um determinado sistema operacional e implementação de JVM? Na verdade, não tenho certeza de qual SO você está usando, porque acabei de verificar minha árvore de desenvolvimento do HotSpot e o Linux usa gettimeofday(&time, NULL)por miliSeconds, o que não é monotônico e não oferece nenhuma garantia de precisão (portanto, em algumas plataformas / kernels exatamente o mesmo problemas como a implementação currentTimeMillis do Windows - para que a pessoa esteja bem ou não esteja). Por outro lado, o nanoTime usa o clock_gettime(CLOCK_MONOTONIC, &tp)que claramente também é a coisa certa a ser usada quando se faz um benchmarking no Linux.
Voo 5/05
Isso nunca causou um problema para mim desde que codifiquei o java em qualquer distro ou kernel do Linux.
Chris Seline
6

Golang versão 45X do python na minha máquina nos códigos Golang abaixo:

package main

import (
"fmt"
"time"
)

const (
n     = 6
iters = 1000
)

var (
x, y, z, w = 34353, 34353, 57768, 1564 //PRNG seeds
)

/* xorshift PRNG
 * Taken from https://en.wikipedia.org/wiki/Xorshift#Example_implementation
 * Used under CC-By-SA */
func myRand() int {
var t uint
t = uint(x ^ (x << 11))
x, y, z = y, z, w
w = int(uint(w^w>>19) ^ t ^ (t >> 8))
return w
}

func main() {
var firstzero, bothzero int
var arr [n + 1]int
var i, j int
x = int(time.Now().Unix())

for i = 0; i < 1<<(n+1); i = i + 1 {
    tmp := i
    for j = 0; j < n+1; j = j + 1 {
        arr[j] = (tmp&1)*(-2) + 1
        tmp >>= 1
    }
    for j = 0; j < iters; j = j + 1 {
        var randArr [n]int
        var flag uint
        var k, first, second int
        for {
            for k = 0; k < n; k = k + 1 {
                randArr[k] = (1 - (myRand() & 3)) % 2
                flag += uint(randArr[k] & 1)
                first += arr[k] * randArr[k]
                second += arr[k+1] * randArr[k]
            }
            if flag != 0 {
                break
            }
        }
        if first == 0 {
            firstzero += 1
            if second == 0 {
                bothzero += 1
            }
        }
    }
}
println("firstzero", firstzero, "bothzero", bothzero)
}

e os códigos python abaixo copiados de cima:

import itertools
import operator
import random

n=6
iters = 1000
firstzero = 0
bothzero = 0

choicesF = filter(any, itertools.product([-1, 0, 0, 1], repeat=n))

for S in itertools.product([-1,1], repeat = n+1):
    for i in xrange(iters):
        F = random.choice(choicesF)
        if not sum(map(operator.mul, F, S[:-1])):
            firstzero += 1
            if not sum(map(operator.mul, F, S[1:])):
                bothzero += 1

print "firstzero", firstzero
print "bothzero", bothzero

e o tempo abaixo:

$time python test.py
firstzero 27349
bothzero 12125

real    0m0.477s
user    0m0.461s
sys 0m0.014s

$time ./hf
firstzero 27253 bothzero 12142

real    0m0.011s
user    0m0.008s
sys 0m0.002s
lunny
fonte
1
você já pensou em usar "github.com/yanatan16/itertools"? Você também diria que isso funcionaria bem em várias goroutines?
ymg
5

C # 0.135s

C # baseado no python simples de Alistair Buxton : 0.278s
C # paralelo: 0.135s
Python da pergunta: 5.907s
python simples de Alistair: 0.853s

Na verdade, não tenho certeza de que essa implementação esteja correta - sua saída é diferente, se você observar os resultados na parte inferior.

Certamente existem algoritmos mais ótimos. Eu apenas decidi usar um algoritmo muito semelhante ao do Python.

C de rosca única

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace ConvolvingArrays
{
    static class Program
    {
        static void Main(string[] args)
        {
            int n=6;
            int iters = 1000;
            int firstzero = 0;
            int bothzero = 0;

            int[] arraySeed = new int[] {-1, 1};
            int[] randomSource = new int[] {-1, 0, 0, 1};
            Random rand = new Random();

            foreach (var S in Enumerable.Repeat(arraySeed, n+1).CartesianProduct())
            {
                for (int i = 0; i < iters; i++)
                {
                    var F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    while (!F.Any(f => f != 0))
                    {
                        F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    }
                    if (Enumerable.Zip(F, S.Take(n), (f, s) => f * s).Sum() == 0)
                    {
                        firstzero++;
                        if (Enumerable.Zip(F, S.Skip(1), (f, s) => f * s).Sum() == 0)
                        {
                            bothzero++;
                        }
                    }
                }
            }

            Console.WriteLine("firstzero {0}", firstzero);
            Console.WriteLine("bothzero {0}", bothzero);
        }

        // itertools.product?
        // http://ericlippert.com/2010/06/28/computing-a-cartesian-product-with-linq/
        static IEnumerable<IEnumerable<T>> CartesianProduct<T>
            (this IEnumerable<IEnumerable<T>> sequences)
        {
            IEnumerable<IEnumerable<T>> emptyProduct =
              new[] { Enumerable.Empty<T>() };
            return sequences.Aggregate(
              emptyProduct,
              (accumulator, sequence) =>
                from accseq in accumulator
                from item in sequence
                select accseq.Concat(new[] { item }));
        }
    }
}

C # paralelo:

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace ConvolvingArrays
{
    static class Program
    {
        static void Main(string[] args)
        {
            int n=6;
            int iters = 1000;
            int firstzero = 0;
            int bothzero = 0;

            int[] arraySeed = new int[] {-1, 1};
            int[] randomSource = new int[] {-1, 0, 0, 1};

            ConcurrentBag<int[]> results = new ConcurrentBag<int[]>();

            // The next line iterates over arrays of length n+1 which contain only -1s and 1s
            Parallel.ForEach(Enumerable.Repeat(arraySeed, n + 1).CartesianProduct(), (S) =>
            {
                int fz = 0;
                int bz = 0;
                ThreadSafeRandom rand = new ThreadSafeRandom();
                for (int i = 0; i < iters; i++)
                {
                    var F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    while (!F.Any(f => f != 0))
                    {
                        F = Enumerable.Range(0, n).Select(_ => randomSource[rand.Next(randomSource.Length)]);
                    }
                    if (Enumerable.Zip(F, S.Take(n), (f, s) => f * s).Sum() == 0)
                    {
                        fz++;
                        if (Enumerable.Zip(F, S.Skip(1), (f, s) => f * s).Sum() == 0)
                        {
                            bz++;
                        }
                    }
                }

                results.Add(new int[] { fz, bz });
            });

            foreach (int[] res in results)
            {
                firstzero += res[0];
                bothzero += res[1];
            }

            Console.WriteLine("firstzero {0}", firstzero);
            Console.WriteLine("bothzero {0}", bothzero);
        }

        // itertools.product?
        // http://ericlippert.com/2010/06/28/computing-a-cartesian-product-with-linq/
        static IEnumerable<IEnumerable<T>> CartesianProduct<T>
            (this IEnumerable<IEnumerable<T>> sequences)
        {
            IEnumerable<IEnumerable<T>> emptyProduct =
              new[] { Enumerable.Empty<T>() };
            return sequences.Aggregate(
              emptyProduct,
              (accumulator, sequence) =>
                from accseq in accumulator
                from item in sequence
                select accseq.Concat(new[] { item }));
        }
    }

    // http://stackoverflow.com/a/11109361/1030702
    public class ThreadSafeRandom
    {
        private static readonly Random _global = new Random();
        [ThreadStatic]
        private static Random _local;

        public ThreadSafeRandom()
        {
            if (_local == null)
            {
                int seed;
                lock (_global)
                {
                    seed = _global.Next();
                }
                _local = new Random(seed);
            }
        }
        public int Next()
        {
            return _local.Next();
        }
        public int Next(int maxValue)
        {
            return _local.Next(maxValue);
        }
    }
}

Saída de teste:

Windows (.NET)

O C # é muito mais rápido no Windows. Provavelmente porque o .NET é mais rápido que o mono.

O tempo do usuário e do sistema não parece funcionar (usado git bashpara o tempo).

$ time /c/Python27/python.exe numpypython.py
firstzero 27413
bothzero 12073

real    0m5.907s
user    0m0.000s
sys     0m0.000s
$ time /c/Python27/python.exe plainpython.py
firstzero 26983
bothzero 12033

real    0m0.853s
user    0m0.000s
sys     0m0.000s
$ time ConvolvingArrays.exe
firstzero 28526
bothzero 6453

real    0m0.278s
user    0m0.000s
sys     0m0.000s
$ time ConvolvingArraysParallel.exe
firstzero 28857
bothzero 6485

real    0m0.135s
user    0m0.000s
sys     0m0.000s

Linux (mono)

bob@phoebe:~/convolvingarrays$ time python program.py
firstzero 27059
bothzero 12131

real    0m11.932s
user    0m11.912s
sys     0m0.012s
bob@phoebe:~/convolvingarrays$ mcs -optimize+ -debug- program.cs
bob@phoebe:~/convolvingarrays$ time mono program.exe
firstzero 28982
bothzero 6512

real    0m1.360s
user    0m1.532s
sys     0m0.872s
bob@phoebe:~/convolvingarrays$ mcs -optimize+ -debug- parallelprogram.cs
bob@phoebe:~/convolvingarrays$ time mono parallelprogram.exe
firstzero 28857
bothzero 6496

real    0m0.851s
user    0m2.708s
sys     0m3.028s
Prumo
fonte
1
Não acho que o código esteja correto como você diz. As saídas não estão corretas.
@Lembik Yea. Eu apreciaria se alguém pudesse me dizer onde está errado - eu não consigo descobrir (apenas ter um entendimento mínimo do que deve fazer não ajuda).
Bob
Seria interessante ver como isso faz com nativo .NET blogs.msdn.com/b/dotnet/archive/2014/04/02/...
Rick Minerich
@Embembik Acabei de analisar tudo isso, tanto quanto posso dizer que deve ser idêntico à outra solução Python ... agora estou realmente confuso.
Bob
4

Haskell: ~ 2000x de aceleração por núcleo

Compile com 'ghc -O3 -funbox-strict-fields -threaded -fllvm' e execute com '+ RTS -Nk' onde k é o número de núcleos em sua máquina.

import Control.Parallel.Strategies
import Data.Bits
import Data.List
import Data.Word
import System.Random

n = 6 :: Int
iters = 1000 :: Int

data G = G !Word !Word !Word !Word deriving (Eq, Show)

gen :: G -> (Word, G)
gen (G x y z w) = let t  = x `xor` (x `shiftL` 11)
                      w' = w `xor` (w `shiftR` 19) `xor` t `xor` (t `shiftR` 8)
                  in (w', G y z w w')  

mask :: Word -> Word
mask = (.&.) $ (2 ^ n) - 1

gen_nonzero :: G -> (Word, G)
gen_nonzero g = let (x, g') = gen g 
                    a = mask x
                in if a == 0 then gen_nonzero g' else (a, g')


data F = F {zeros  :: !Word, 
            posneg :: !Word} deriving (Eq, Show)

gen_f :: G -> (F, G)       
gen_f g = let (a, g')  = gen_nonzero g
              (b, g'') = gen g'
          in  (F a $ mask b, g'')

inner :: Word -> F -> Int
inner s (F zs pn) = let s' = complement $ s `xor` pn
                        ones = s' .&. zs
                        negs = (complement s') .&. zs
                    in popCount ones - popCount negs

specialised_convolve :: Word -> F -> (Int, Int)
specialised_convolve s f@(F zs pn) = (inner s f', inner s f) 
    where f' = F (zs `shiftL` 1) (pn `shiftL` 1)

ss :: [Word]
ss = [0..2 ^ (n + 1) - 1]

main_loop :: [G] -> (Int, Int)
main_loop gs = foldl1' (\(fz, bz) (fz', bz') -> (fz + fz', bz + bz')) . parMap rdeepseq helper $ zip ss gs
    where helper (s, g) = go 0 (0, 0) g
                where go k u@(fz, bz) g = if k == iters 
                                              then u 
                                              else let (f, g') = gen_f g
                                                       v = case specialised_convolve s f
                                                               of (0, 0) -> (fz + 1, bz + 1)
                                                                  (0, _) -> (fz + 1, bz)
                                                                  _      -> (fz, bz)
                                                   in go (k + 1) v g'

seed :: IO G                                        
seed = do std_g <- newStdGen
          let [x, y, z, w] = map fromIntegral $ take 4 (randoms std_g :: [Int])
          return $ G x y z w

main :: IO ()
main = (sequence $ map (const seed) ss) >>= print . main_loop
user1502040
fonte
2
Então, com 4 núcleos , são mais de 9000 ?! Não há como isso estar certo.
Cees Timmerman
A lei de Amdahl declara que a aceleração da paralelização não é linear ao número de unidades de processamento paralelo. em vez disso, eles apenas fornecem retornos
escuros
@xaedes A aceleração parece essencialmente linear para um número baixo de núcleos
user1502040 12/07
3

Rubi

Ruby (2.1.0) 0.277s
Ruby (2.1.1) 0.281s
Python (Alistair Buxton) 0.330s
Python (alemão) 0.097s

n = 6
iters = 1000
first_zero = 0
both_zero = 0

choices = [-1, 0, 0, 1].repeated_permutation(n).select{|v| [0] != v.uniq}

def convolve(v1, v2)
  [0, 1].map do |i|
    r = 0
    6.times do |j|
      r += v1[i+j] * v2[j]
    end
    r
  end
end

[-1, 1].repeated_permutation(n+1) do |s|
  iters.times do
    f = choices.sample
    fs = convolve s, f
    if 0 == fs[0]
      first_zero += 1
      if 0 == fs[1]
        both_zero += 1
      end
    end
  end
end

puts 'firstzero %i' % first_zero
puts 'bothzero %i' % both_zero
Landstander
fonte
3

thread não estaria completo sem PHP

6.6x mais rápido

PHP v5.5.9 - 1.223 0,646 seg;

vs

Python v2.7.6 - 8.072 sec

<?php

$n = 6;
$iters = 1000;
$firstzero = 0;
$bothzero = 0;

$x=time();
$y=34353;
$z=57768;
$w=1564; //PRNG seeds

function myRand() {
    global $x;
    global $y;
    global $z;
    global $w;
    $t = $x ^ ($x << 11);
    $x = $y; $y = $z; $z = $w;
    return $w = $w ^ ($w >> 19) ^ $t ^ ($t >> 8);
}

function array_cartesian() {
    $_ = func_get_args();
    if (count($_) == 0)
        return array();
    $a = array_shift($_);
    if (count($_) == 0)
        $c = array(array());
    else
        $c = call_user_func_array(__FUNCTION__, $_);
    $r = array();
    foreach($a as $v)
        foreach($c as $p)
            $r[] = array_merge(array($v), $p);
    return $r;
}

function rand_array($a, $n)
{
    $r = array();
    for($i = 0; $i < $n; $i++)
        $r[] = $a[myRand()%count($a)];
    return $r;
}

function convolve($a, $b)
{
    // slows down
    /*if(count($a) < count($b))
        return convolve($b,$a);*/
    $result = array();
    $w = count($a) - count($b) + 1;
    for($i = 0; $i < $w; $i++){
        $r = 0;
        for($k = 0; $k < count($b); $k++)
            $r += $b[$k] * $a[$i + $k];
        $result[] = $r;
    }
    return $result;
}

$cross = call_user_func_array('array_cartesian',array_fill(0,$n+1,array(-1,1)));

foreach($cross as $S)
    for($i = 0; $i < $iters; $i++){
        while(true)
        {
            $F = rand_array(array(-1,0,0,1), $n);
            if(in_array(-1, $F) || in_array(1, $F))
                break;
        }
        $FS = convolve($S, $F);
        if(0==$FS[0]) $firstzero += 1;
        if(0==$FS[0] && 0==$FS[1]) $bothzero += 1;
    }

echo "firstzero $firstzero\n";
echo "bothzero $bothzero\n";
  • Usava um gerador aleatório personalizado (roubado da resposta C), o PHP um é péssimo e os números não correspondem
  • convolve função simplificada um pouco para ser mais rápido
  • A verificação de array apenas com zeros também é muito otimizada (consulte $Fe $FSverificações).

Saídas:

$ time python num.py 
firstzero 27050
bothzero 11990

real    0m8.072s
user    0m8.037s
sys 0m0.024s
$ time php num.php
firstzero 27407
bothzero 12216

real    0m1.223s
user    0m1.210s
sys 0m0.012s

Editar. A segunda versão do script funciona apenas para 0.646 sec:

<?php

$n = 6;
$iters = 1000;
$firstzero = 0;
$bothzero = 0;

$x=time();
$y=34353;
$z=57768;
$w=1564; //PRNG seeds

function myRand() {
    global $x;
    global $y;
    global $z;
    global $w;
    $t = $x ^ ($x << 11);
    $x = $y; $y = $z; $z = $w;
    return $w = $w ^ ($w >> 19) ^ $t ^ ($t >> 8);
}

function array_cartesian() {
    $_ = func_get_args();
    if (count($_) == 0)
        return array();
    $a = array_shift($_);
    if (count($_) == 0)
        $c = array(array());
    else
        $c = call_user_func_array(__FUNCTION__, $_);
    $r = array();
    foreach($a as $v)
        foreach($c as $p)
            $r[] = array_merge(array($v), $p);
    return $r;
}

function convolve($a, $b)
{
    // slows down
    /*if(count($a) < count($b))
        return convolve($b,$a);*/
    $result = array();
    $w = count($a) - count($b) + 1;
    for($i = 0; $i < $w; $i++){
        $r = 0;
        for($k = 0; $k < count($b); $k++)
            $r += $b[$k] * $a[$i + $k];
        $result[] = $r;
    }
    return $result;
}

$cross = call_user_func_array('array_cartesian',array_fill(0,$n+1,array(-1,1)));

$choices = call_user_func_array('array_cartesian',array_fill(0,$n,array(-1,0,0,1)));

foreach($cross as $S)
    for($i = 0; $i < $iters; $i++){
        while(true)
        {
            $F = $choices[myRand()%count($choices)];
            if(in_array(-1, $F) || in_array(1, $F))
                break;
        }
        $FS = convolve($S, $F);
        if(0==$FS[0]){
            $firstzero += 1;
            if(0==$FS[1])
                $bothzero += 1;
        }
    }

echo "firstzero $firstzero\n";
echo "bothzero $bothzero\n";
Vitaly Dyatlov
fonte
3

Solução F #

O tempo de execução é de 0,030s quando compilado para x86 no CLR Core i7 4 (8) a 3,4 Ghz

Não faço ideia se o código está correto.

  • Otimização funcional (dobra em linha) -> 0,026s
  • Construindo via projeto do console -> 0.022s
  • Adicionado um algoritmo melhor para a geração de matrizes de permutação -> 0,018s
  • Mono para Windows -> 0.089s
  • Executando o script Python de Alistair -> 0.259s
let inline ffoldi n f state =
    let mutable state = state
    for i = 0 to n - 1 do
        state <- f state i
    state

let product values n =
    let p = Array.length values
    Array.init (pown p n) (fun i ->
        (Array.zeroCreate n, i)
        |> ffoldi n (fun (result, i') j ->
            result.[j] <- values.[i' % p]
            result, i' / p
        )
        |> fst
    )

let convolute signals filter =
    let m = Array.length signals
    let n = Array.length filter
    let len = max m n - min m n + 1

    Array.init len (fun offset ->
        ffoldi n (fun acc i ->
            acc + filter.[i] * signals.[m - 1 - offset - i]
        ) 0
    )

let n = 6
let iters = 1000

let next =
    let arrays =
        product [|-1; 0; 0; 1|] n
        |> Array.filter (Array.forall ((=) 0) >> not)
    let rnd = System.Random()
    fun () -> arrays.[rnd.Next arrays.Length]

let signals = product [|-1; 1|] (n + 1)

let firstzero, bothzero =
    ffoldi signals.Length (fun (firstzero, bothzero) i ->
        let s = signals.[i]
        ffoldi iters (fun (first, both) _ ->
            let f = next()
            match convolute s f with
            | [|0; 0|] -> first + 1, both + 1
            | [|0; _|] -> first + 1, both
            | _ -> first, both
        ) (firstzero, bothzero)
    ) (0, 0)

printfn "firstzero %i" firstzero
printfn "bothzero %i" bothzero
David Grenier
fonte
2

Q, 0,296 seg

n:6; iter:1000  /parametrization (constants)
c:n#0           /auxiliar constant (sequence 0 0.. 0 (n))
A:B:();         /A and B accumulates results of inner product (firstresult, secondresult)

/S=sequence with all arrays of length n+1 with values -1 and 1
S:+(2**m)#/:{,/x#/:-1 1}'m:|n(2*)\1 

f:{do[iter; F:c; while[F~c; F:n?-1 0 0 1]; A,:+/F*-1_x; B,:+/F*1_x];} /hard work
f'S               /map(S,f)
N:~A; +/'(N;N&~B) / ~A is not A (or A=0) ->bitmap.  +/ is sum (population over a bitmap)
                  / +/'(N;N&~B) = count firstResult=0, count firstResult=0 and secondResult=0

Q é uma linguagem orientada a coleções (kx.com)

Código reescrito para explorar Q idiomático, mas nenhuma outra otimização inteligente

Linguagens de script otimizam o tempo do programador, não o tempo de execução

  • Q não é a melhor ferramenta para este problema

Primeira tentativa de codificação = não é um vencedor, mas um tempo razoável (aprox. 30x de aceleração)

  • bastante competitivo entre intérpretes
  • pare e escolha outro problema

NOTAS.-

  • programa usa semente padrão (execs repetíveis) Para escolher outra semente para uso aleatório de gerador \S seed
  • O resultado é dado como um squence de duas polegadas, portanto, existe um sufixo i final no segundo valor 27421 12133i -> lido como (27241, 12133)
  • Tempo sem contar a inicialização do intérprete. \t sentence mede o tempo consumido por essa frase
J. Sendra
fonte
Muito interessante obrigado.
1

Julia: 12.149 6.929 s

Apesar das reivindicações de velocidade , o tempo inicial de compilação do JIT nos impede!

Observe que o código Julia a seguir é efetivamente uma tradução direta do código Python original (sem otimizações) como uma demonstração de que você pode transferir facilmente sua experiência de programação para uma linguagem mais rápida;)

require("Iterators")

n = 6
iters = 1000
firstzero = 0
bothzero = 0

for S in Iterators.product(fill([-1,1], n+1)...)
    for i = 1:iters
        F = [[-1 0 0 1][rand(1:4)] for _ = 1:n]
        while all((x) -> round(x,8) == 0, F)
            F = [[-1 0 0 1][rand(1:4)] for _ = 1:n]
        end
        FS = conv(F, [S...])
        if round(FS[1],8) == 0
            firstzero += 1
        end
        if all((x) -> round(x,8) == 0, FS)
            bothzero += 1
        end
    end
end

println("firstzero ", firstzero)
println("bothzero ", bothzero)

Editar

Correr com n = 8leva 32.935 s. Considerando que a complexidade desse algoritmo é O(2^n), então 4 * (12.149 - C) = (32.935 - C), Cuma constante que representa o tempo de compilação do JIT. Resolvendo Cisso C = 5.2203, sugerimos que o tempo real de execução n = 6é de 6.929 s.

ágar ágil
fonte
Que tal aumentar n para 8 para ver se Julia se destaca então?
Isso ignora muitas dicas de desempenho aqui: julia.readthedocs.org/en/latest/manual/performance-tips . Veja também a outra entrada de Julia, que se sai significativamente melhor. A submissão é apreciado embora :-)
StefanKarpinski
0

Ferrugem, 6,6 ms, aceleração de 1950x

Praticamente uma tradução direta do código de Alistair Buxton para Rust. Pensei em usar vários núcleos com rayon (concorrência sem medo!), Mas isso não melhorou o desempenho, provavelmente porque já é muito rápido.

extern crate itertools;
extern crate rand;
extern crate time;

use itertools::Itertools;
use rand::{prelude::*, prng::XorShiftRng};
use std::iter;
use time::precise_time_ns;

fn main() {
    let start = precise_time_ns();

    let n = 6;
    let iters = 1000;
    let mut first_zero = 0;
    let mut both_zero = 0;
    let choices_f: Vec<Vec<i8>> = iter::repeat([-1, 0, 0, 1].iter().cloned())
        .take(n)
        .multi_cartesian_product()
        .filter(|i| i.iter().any(|&x| x != 0))
        .collect();
    // xorshift RNG is faster than default algorithm designed for security
    // rather than performance.
    let mut rng = XorShiftRng::from_entropy(); 
    for s in iter::repeat(&[-1, 1]).take(n + 1).multi_cartesian_product() {
        for _ in 0..iters {
            let f = rng.choose(&choices_f).unwrap();
            if f.iter()
                .zip(&s[..s.len() - 1])
                .map(|(a, &b)| a * b)
                .sum::<i8>() == 0
            {
                first_zero += 1;
                if f.iter().zip(&s[1..]).map(|(a, &b)| a * b).sum::<i8>() == 0 {
                    both_zero += 1;
                }
            }
        }
    }
    println!("first_zero = {}\nboth_zero = {}", first_zero, both_zero);

    println!("runtime {} ns", precise_time_ns() - start);
}

E Cargo.toml, como eu uso dependências externas:

[package]
name = "how_slow_is_python"
version = "0.1.0"

[dependencies]
itertools = "0.7.8"
rand = "0.5.3"
time = "0.1.40"

Comparação de velocidade:

$ time python2 py.py
firstzero: 27478
bothzero: 12246
12.80user 0.02system 0:12.90elapsed 99%CPU (0avgtext+0avgdata 23328maxresident)k
0inputs+0outputs (0major+3544minor)pagefaults 0swaps
$ time target/release/how_slow_is_python
first_zero = 27359
both_zero = 12162
runtime 6625608 ns
0.00user 0.00system 0:00.00elapsed 100%CPU (0avgtext+0avgdata 2784maxresident)k
0inputs+0outputs (0major+189minor)pagefaults 0swaps

6625608 ns é de cerca de 6,6 ms. Isso significa aceleração de 1950 vezes. Existem muitas otimizações possíveis aqui, mas eu estava buscando mais legibilidade do que desempenho. Uma otimização possível seria usar matrizes em vez de vetores para armazenar opções, pois elas sempre terão nelementos. Também é possível usar o RNG que não seja o XorShift, pois, embora o Xorshift seja mais rápido que o HC-128 CSPRNG padrão, é mais lento que o ingênuo dos algoritmos PRNG.

Konrad Borowski
fonte