Encontrar correlaciones aproximadas

14

Considere una cadena binaria Sde longitud n. Indexando desde 1, podemos calcular las distancias de Hamming entre S[1..i+1]y S[n-i..n]para todos ien orden de 0a n-1. La distancia de Hamming entre dos cadenas de igual longitud es el número de posiciones en las que los símbolos correspondientes son diferentes. Por ejemplo,

S = 01010

da

[0, 2, 0, 4, 0].

Esto se debe a que 0coincide 0, 01tiene una distancia de Hamming de dos a 10, 010coincide 010, 0101tiene una distancia de Hamming de cuatro a 1010 y finalmente 01010coincide.

Sin embargo, solo nos interesan las salidas donde la distancia de Hamming es como máximo 1. Entonces, en esta tarea, informaremos Ysi la distancia de Hamming es como máximo una y Notra. Entonces, en nuestro ejemplo anterior, obtendríamos

[Y, N, Y, N, Y]

Definir f(n)a ser el número de matrices distintas de Ys y Ns se obtiene cuando se itera sobre todas las 2^ndiferentes cadenas de bits posibles Sde longitud n.

Tarea

Para aumentar a npartir de 1, su código debería salir f(n).

Ejemplo de respuestas

Para n = 1..24, las respuestas correctas son:

1, 1, 2, 4, 6, 8, 14, 18, 27, 36, 52, 65, 93, 113, 150, 188, 241, 279, 377, 427, 540, 632, 768, 870

Puntuación

Su código debe iterar desde n = 1dar la respuesta para cada uno npor turno. Voy a cronometrar toda la carrera, matándola después de dos minutos.

Tu puntaje es el más alto n que alcanzas en ese momento.

En caso de empate, la primera respuesta gana.

¿Dónde se probará mi código?

Ejecutaré su código en mi computadora portátil Windows 7 (un poco vieja) con cygwin. Como resultado, brinde toda la ayuda que pueda para ayudar a que esto sea más fácil.

Mi laptop tiene 8GB de RAM y una CPU Intel i7 [email protected] GHz (Broadwell) con 2 núcleos y 4 hilos. El conjunto de instrucciones incluye SSE4.2, AVX, AVX2, FMA3 y TSX.

Entradas principales por idioma

  • n = 40 en Rust usando CryptoMiniSat, por Anders Kaseorg. (En la máquina virtual invitada de Lubuntu en Vbox).
  • n = 35 en C ++ usando la biblioteca BuDDy, por Christian Seviers. (En la máquina virtual invitada de Lubuntu en Vbox).
  • n = 34 en Clingo por Anders Kaseorg. (En la máquina virtual invitada de Lubuntu en Vbox).
  • n = 31 en Rust por Anders Kaseorg.
  • n = 29 en Clojure por NikoNyrh.
  • n = 29 en C por bartavelle.
  • n = 27 en Haskell por bartavelle
  • n = 24 en Pari / gp por alephalpha.
  • n = 22 en Python 2 + pypy por mí.
  • n = 21 en Mathematica por alephalpha. (Autoinformado)

Recompensas futuras

Ahora daré una recompensa de 200 puntos por cualquier respuesta que llegue a n = 80 en mi máquina en dos minutos.


fuente
¿Conoces algún truco que permita a alguien encontrar un algoritmo más rápido que una fuerza bruta ingenua? Si no, este desafío es "implemente esto en x86" (o tal vez si conocemos su GPU ...).
Jonathan Allan
@JonathanAllan Ciertamente es posible acelerar un enfoque muy ingenuo. No sé exactamente qué tan rápido puede llegar. Curiosamente, si cambiamos la pregunta para que obtenga una S si la distancia de Hamming es como máximo 0 y una N de lo contrario, entonces hay una fórmula de forma cerrada conocida.
@Lembik ¿Medimos el tiempo de CPU o el tiempo real?
flawr
@flawr Estoy midiendo el tiempo real, pero lo ejecuté varias veces y tomé el mínimo para eliminar las rarezas.

Respuestas:

9

Rust + CryptoMiniSat , n ° 41

src/main.rs

extern crate cryptominisat;
extern crate itertools;

use std::iter::once;
use cryptominisat::{Lbool, Lit, Solver};
use itertools::Itertools;

fn make_solver(n: usize) -> (Solver, Vec<Lit>) {
    let mut solver = Solver::new();
    let s: Vec<Lit> = (1..n).map(|_| solver.new_var()).collect();
    let d: Vec<Vec<Lit>> = (1..n - 1)
        .map(|k| {
                 (0..n - k)
                     .map(|i| (if i == 0 { s[k - 1] } else { solver.new_var() }))
                     .collect()
             })
        .collect();
    let a: Vec<Lit> = (1..n - 1).map(|_| solver.new_var()).collect();
    for k in 1..n - 1 {
        for i in 1..n - k {
            solver.add_xor_literal_clause(&[s[i - 1], s[k + i - 1], d[k - 1][i]], true);
        }
        for t in (0..n - k).combinations(2) {
            solver.add_clause(&t.iter()
                                   .map(|&i| d[k - 1][i])
                                   .chain(once(!a[k - 1]))
                                   .collect::<Vec<_>>()
                                   [..]);
        }
        for t in (0..n - k).combinations(n - k - 1) {
            solver.add_clause(&t.iter()
                                   .map(|&i| !d[k - 1][i])
                                   .chain(once(a[k - 1]))
                                   .collect::<Vec<_>>()
                                   [..]);
        }
    }
    (solver, a)
}

fn search(n: usize,
          solver: &mut Solver,
          a: &Vec<Lit>,
          assumptions: &mut Vec<Lit>,
          k: usize)
          -> usize {
    match solver.solve_with_assumptions(assumptions) {
        Lbool::True => search_sat(n, solver, a, assumptions, k),
        Lbool::False => 0,
        Lbool::Undef => panic!(),
    }
}

fn search_sat(n: usize,
              solver: &mut Solver,
              a: &Vec<Lit>,
              assumptions: &mut Vec<Lit>,
              k: usize)
              -> usize {
    if k >= n - 1 {
        1
    } else {
        let s = solver.is_true(a[k - 1]);
        assumptions.push(if s { a[k - 1] } else { !a[k - 1] });
        let c = search_sat(n, solver, a, assumptions, k + 1);
        assumptions.pop();
        assumptions.push(if s { !a[k - 1] } else { a[k - 1] });
        let c1 = search(n, solver, a, assumptions, k + 1);
        assumptions.pop();
        c + c1
    }
}

fn f(n: usize) -> usize {
    let (mut solver, proj) = make_solver(n);
    search(n, &mut solver, &proj, &mut vec![], 1)
}

fn main() {
    for n in 1.. {
        println!("{}: {}", n, f(n));
    }
}

Cargo.toml

[package]
name = "correlations-cms"
version = "0.1.0"
authors = ["Anders Kaseorg <[email protected]>"]

[dependencies]
cryptominisat = "5.0.1"
itertools = "0.6.0"

Cómo funciona

Esto hace una búsqueda recursiva a través del árbol de todas las asignaciones parciales a los prefijos de la matriz Y / N, utilizando un solucionador SAT para verificar en cada paso si la asignación parcial actual es consistente y podar la búsqueda si no. CryptoMiniSat es el solucionador SAT adecuado para este trabajo debido a sus optimizaciones especiales para las cláusulas XOR.

Las tres familias de restricciones son

S iS k + iD ki , para 1 ≤ kn - 2, 0 ≤ i ≤ n - k ;
D ki 1D ki 2 ∨ ¬ A k , para 1 ≤ kn - 2, 0 ≤ i 1 < i 2n - k ;
¬ D ki 1 ∨ ⋯ ∨ ¬ D ki n - k - 1A k , para 1 ≤ kn - 2, 0 ≤ i 1 <⋯ < i n - k - 1n - k ;

excepto que, como optimización, S 0 se fuerza a falso, de modo que D k 0 es simplemente igual a S k .

Anders Kaseorg
fuente
2
Woohoooooo! :)
Todavía estoy tratando de compilar esto en Windows (usando cygwin + gcc). Cloné cryptominisat y lo compilé. Pero todavía no sé cómo compilar el código de óxido. Cuando lo hago cargo build, recibo--- stderr CMake Error: Could not create named generator Visual Studio 14 2015 Win64
2
@ rahnema1 Gracias, pero parece que el problema es con el sistema de compilación CMake de la biblioteca C ++ incrustada en la caja de cryptominisat, no con Rust en sí.
Anders Kaseorg
1
@Lembik, obtengo un 404 de esa pasta.
Mego
1
@ChristianSievers Buena pregunta. Eso funciona, pero parece ser un poco más lento (2 × más o menos). No estoy seguro de por qué no debería ser tan bueno, por lo que tal vez CryptoMiniSat no haya sido optimizado para ese tipo de carga de trabajo incremental.
Anders Kaseorg
9

Óxido, n ≈ 30 o 31 o 32

En mi computadora portátil (dos núcleos, i5-6200U), esto pasa por n = 1, ..., 31 en 53 segundos, usando aproximadamente 2.5 GiB de memoria, o por n = 1, ..., 32 en 105 segundos, usando aproximadamente 5 GiB de la memoria Compilar cargo build --releasey ejecutar target/release/correlations.

src/main.rs

extern crate rayon;

type S = u32;
const S_BITS: u32 = 32;

fn cat(mut a: Vec<S>, mut b: Vec<S>) -> Vec<S> {
    if a.capacity() >= b.capacity() {
        a.append(&mut b);
        a
    } else {
        b.append(&mut a);
        b
    }
}

fn search(n: u32, i: u32, ss: Vec<S>) -> u32 {
    if ss.is_empty() {
        0
    } else if 2 * i + 1 > n {
        search_end(n, i, ss)
    } else if 2 * i + 1 == n {
        search2(n, i, ss.into_iter().flat_map(|s| vec![s, s | 1 << i]))
    } else {
        search2(n,
                i,
                ss.into_iter()
                    .flat_map(|s| {
                                  vec![s,
                                       s | 1 << i,
                                       s | 1 << n - i - 1,
                                       s | 1 << i | 1 << n - i - 1]
                              }))
    }
}

fn search2<SS: Iterator<Item = S>>(n: u32, i: u32, ss: SS) -> u32 {
    let (shift, mask) = (n - i - 1, !(!(0 as S) << i + 1));
    let close = |s: S| {
        let x = (s ^ s >> shift) & mask;
        x & x.wrapping_sub(1) == 0
    };
    let (ssy, ssn) = ss.partition(|&s| close(s));
    let (cy, cn) = rayon::join(|| search(n, i + 1, ssy), || search(n, i + 1, ssn));
    cy + cn
}

fn search_end(n: u32, i: u32, ss: Vec<S>) -> u32 {
    if i >= n - 1 { 1 } else { search_end2(n, i, ss) }
}

fn search_end2(n: u32, i: u32, mut ss: Vec<S>) -> u32 {
    let (shift, mask) = (n - i - 1, !(!(0 as S) << i + 1));
    let close = |s: S| {
        let x = (s ^ s >> shift) & mask;
        x & x.wrapping_sub(1) == 0
    };
    match ss.iter().position(|&s| close(s)) {
        Some(0) => {
            match ss.iter().position(|&s| !close(s)) {
                Some(p) => {
                    let (ssy, ssn) = ss.drain(p..).partition(|&s| close(s));
                    let (cy, cn) = rayon::join(|| search_end(n, i + 1, cat(ss, ssy)),
                                               || search_end(n, i + 1, ssn));
                    cy + cn
                }
                None => search_end(n, i + 1, ss),
            }
        }
        Some(p) => {
            let (ssy, ssn) = ss.drain(p..).partition(|&s| close(s));
            let (cy, cn) = rayon::join(|| search_end(n, i + 1, ssy),
                                       || search_end(n, i + 1, cat(ss, ssn)));
            cy + cn
        }
        None => search_end(n, i + 1, ss),
    }
}

fn main() {
    for n in 1..S_BITS + 1 {
        println!("{}: {}", n, search(n, 1, vec![0, 1]));
    }
}

Cargo.toml

[package]
name = "correlations"
version = "0.1.0"
authors = ["Anders Kaseorg <[email protected]>"]

[dependencies]
rayon = "0.7.0"

Pruébalo en línea!

También tengo una variante un poco más lenta que usa mucha menos memoria.

Anders Kaseorg
fuente
¿Qué optimizaciones has usado?
1
@Lembik La mayor optimización, además de hacer todo con aritmética bit a bit en un lenguaje compilado, es usar solo la cantidad de no determinismo que sea necesario para definir un prefijo de la matriz Y / N. Hago una búsqueda recursiva en posibles prefijos de la matriz Y / N, tomando un vector de posibles cadenas que alcanzan ese prefijo, pero solo las cadenas cuyo centro no examinado está lleno de ceros. Dicho esto, esta sigue siendo una búsqueda exponencial, y estas optimizaciones solo lo aceleran por factores polinómicos.
Anders Kaseorg
Es una buena respuesta. Gracias. Espero que alguien profundice en la combinatoria para acelerar significativamente.
@Lembik He arreglado un error de pérdida de memoria, he hecho más micro-optimización y he agregado paralelismo. Vuelva a realizar la prueba cuando tenga la oportunidad. Espero aumentar mi puntaje en 1 o 2. ¿Tiene ideas combinatorias en mente para mayores aceleraciones? No se me ocurrió nada.
Anders Kaseorg
1
@Lembik No hay una fórmula dada en la entrada OEIS. (El código de Mathematica allí también parece usar fuerza bruta). Si conoce uno, es posible que desee contarles al respecto.
Christian Sievers
6

C ++ usando la biblioteca BuDDy

Un enfoque diferente: tener una fórmula binaria (como diagrama de decisión binario ) que toma los bits de Scomo entrada y es cierto si da algunos valores fijos de Yo Nen ciertas posiciones seleccionadas. Si esa fórmula no es constante falsa, seleccione una posición libre y recurse, probando ambas YyN . Si no hay una posición libre, hemos encontrado un posible valor de salida. Si la fórmula es constante falsa, retroceda.

Esto funciona relativamente razonable porque hay muy pocos valores posibles, por lo que a menudo podemos retroceder temprano. Intenté una idea similar con un solucionador SAT, pero fue menos exitosa.

#include<vector>
#include<iostream>
#include<bdd.h>

// does vars[0..i-1] differ from vars[n-i..n-1] in at least two positions?
bdd cond(int i, int n, const std::vector<bdd>& vars){
  bdd x1 { bddfalse };
  bdd xs { bddfalse };
  for(int k=0; k<i; ++k){
    bdd d { vars[k] ^ vars[n-i+k] };
    xs |= d & x1;
    x1 |= d;
  }
  return xs;
}

void expand(int i, int n, int &c, const std::vector<bdd>& conds, bdd x){
  if (x==bddfalse)
    return;
  if (i==n-2){
    ++c;
    return;
  }

  expand(i+1,n,c,conds, x & conds[2*i]);
  x &= conds[2*i+1];
  expand(i+1,n,c,conds, x);
}

int count(int n){
  if (n==1)   // handle trivial case
    return 1;
  bdd_setvarnum(n-1);
  std::vector<bdd> vars {};
  vars.push_back(bddtrue); // assume first bit is 1
  for(int i=0; i<n-1; ++i)
    if (i%2==0)            // vars in mixed order
      vars.push_back(bdd_ithvar(i/2));
    else
      vars.push_back(bdd_ithvar(n-2-i/2));
  std::vector<bdd> conds {};
  for(int i=n-1; i>1; --i){ // handle long blocks first
    bdd cnd { cond(i,n,vars) };
    conds.push_back( cnd );
    conds.push_back( !cnd );
  }
  int c=0;
  expand(0,n,c,conds,bddtrue);
  return c;
}

int main(void){
  bdd_init(20000000,1000000);
  bdd_gbc_hook(nullptr); // comment out to see GC messages
  for(int n=1; ; ++n){
    std::cout << n << " " << count(n) << "\n" ;
  }
}

Para compilar con debian 8 (jessie), instálalo libbdd-devy hazlo g++ -std=c++11 -O3 -o hb hb.cpp -lbdd. Puede ser útil aumentar el primer argumento a bdd_initmás.

Christian Sievers
fuente
Esto se ve interesante. ¿A qué llegas con esto?
@Lembik Tengo 31 en 100 en hardware muy antiguo que no me deja responder más rápido
Christian Sievers
Cualquier ayuda que pueda dar sobre cómo compilar esto en Windows (por ejemplo, usando cygwin) recibió con gratitud.
@Lembik No sé sobre Windws pero github.com/fd00/yacp/tree/master/buddy parece útil wrt cygwin
Christian Sievers
1
Wow, está bien, me has convencido de que necesito agregar esta biblioteca a mi kit de herramientas. ¡Bien hecho!
Anders Kaseorg
4

Clingo, n. ° 30 o 31 34

Me sorprendió un poco ver que cinco líneas de código de Clingo superan mi solución Rust de fuerza bruta y se acercan mucho a la solución BuDDy de Christian: parece que también superaría eso con un límite de tiempo más alto.

corr.lp

{s(2..n)}.
d(K,I) :- K=1..n-2, I=1..n-K, s(I), not s(K+I).
d(K,I) :- K=1..n-2, I=1..n-K, not s(I), s(K+I).
a(K) :- K=1..n-2, {d(K,1..n-K)} 1.
#show a/1.

corr.sh

#!/bin/bash
for ((n=1;;n++)); do
    echo "$n $(clingo corr.lp --project -q -n 0 -c n=$n | sed -n 's/Models *: //p')"
done

trama

Anders Kaseorg
fuente
¡Esto es genial! Parece de su gráfico que la solución BuDDy empeora repentinamente. ¿Alguna idea de por qué?
@Lembik No he investigado BuDDy lo suficiente como para estar seguro, pero ¿tal vez se quede sin memoria caché en ese momento?
Anders Kaseorg
¡Guauu! Creo que un primer valor más alto bdd_initpodría ayudar, o permitir aumentar más la tabla de nodos al llamar bdd_setmaxincreasecon un valor muy por encima del valor predeterminado de 50000. - ¿Está utilizando la versión modificada de mi programa?
Christian Sievers
2
Me encanta tu gráfico.
1
Obtienes un sorprendente aumento del rendimiento al usar la opción --configuration=crafty( jumpyy trendydar resultados similares).
Christian Sievers
2

Pari / GP , 23

Por defecto, Pari / GP limita su tamaño de pila a 8 MB. La primera línea del código default(parisize, "4g"), establece este límite en 4 GB. Si todavía da un flujo de stackoverflow, puede configurarlo en 8 GB.

default(parisize, "4g")
f(n) = #vecsort([[2 > hammingweight(bitxor(s >> (n-i) , s % 2^i)) | i <- [2..n-1]] | s <- [0..2^(n-1)]], , 8)
for(n = 1, 100, print(n " -> " f(n)))
alephalpha
fuente
Alcanza 22 y luego da un stackoverflow.
Llega a 24 ahora.
2

Clojure, 29 en 75 38 segundos, 30 en 80 y 31 en 165

Tiempos de ejecución de Intel i7 6700K , el uso de memoria es inferior a 200 MB.

project.clj (usa com.climate / claypoole para multihilo):

(defproject tests "0.0.1-SNAPSHOT"
  :description "FIXME: write description"
  :dependencies [[org.clojure/clojure "1.8.0"]
                 [com.climate/claypoole "1.1.4"]]
  :javac-options ["-target" "1.6" "-source" "1.6" "-Xlint:-options"]
  :aot [tests.core]
  :main tests.core)

Código fuente:

(ns tests.core
  (:require [com.climate.claypoole :as cp]
            [clojure.set])
  (:gen-class))

(defn f [N]
  (let [n-threads   (.. Runtime getRuntime availableProcessors)
        mask-offset (- 31 N)
        half-N      (quot N 2)
        mid-idx     (bit-shift-left 1 half-N)
        end-idx     (bit-shift-left 1 (dec N))
        lower-half  (bit-shift-right 0x7FFFFFFF mask-offset)
        step        (bit-shift-left 1 12)
        bitcount
          (fn [n]
            (loop [i 0 result 0]
              (if (= i N)
                result
                (recur
                  (inc i)
                  (-> n
                      (bit-xor (bit-shift-right n i))
                      (bit-and (bit-shift-right 0x7FFFFFFF (+ mask-offset i)))
                      Integer/bitCount
                      (< 2)
                      (if (+ result (bit-shift-left 1 i))
                          result))))))]
    (->>
      (cp/upfor n-threads [start (range 0 end-idx step)]
        (->> (for [i      (range start (min (+ start step) end-idx))
                   :when  (<= (Integer/bitCount (bit-shift-right i mid-idx))
                              (Integer/bitCount (bit-and         i lower-half)))]
               (bitcount i))
             (into #{})))
      (reduce clojure.set/union)
      count)))

(defn -main [n]
  (let [n-iters 5]
    (println "Calculating f(n) from 1 to" n "(inclusive)" n-iters "times")
    (doseq [i (range n-iters)]
      (->> n read-string inc (range 1) (map f) doall println time)))
  (shutdown-agents)
  (System/exit 0))

Una solución de fuerza bruta, cada subproceso pasa por un subconjunto del rango (2 ^ 12 elementos) y crea un conjunto de valores enteros que indican patrones detectados. Estos se "unen" juntos y, por lo tanto, se calcula el recuento distinto. Espero que el código no sea demasiado complicado de seguir aunque use muchas macros de subprocesos . Mimain ejecuta la prueba varias veces para calentar JVM.

Actualización: Iterando solo la mitad de los enteros, obtiene el mismo resultado debido a la simetría. También se omiten los números con un mayor número de bits en la mitad inferior del número, ya que también producen duplicados.

Uberjar preconstruido ( v1 ) (3.7 MB):

$ wget https://s3-eu-west-1.amazonaws.com/nikonyrh-public/misc/so-124424-v2.jar
$ java -jar so-124424-v2.jar 29
Calculating f(n) from 1 to 29 (inclusive) 5 times
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 41341.863703 msecs"
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 37752.118265 msecs"
(1 1 2 4 6 8 14 18 27 36 52 65 93 113 150 188 241 279 377 427 540 632 768 870 1082 1210 1455 1656 1974)
"Elapsed time: 38568.406528 msecs"
[ctrl+c]

Resultados en diferentes hardwares, tiempo de ejecución esperado es O(n * 2^n)?

i7-6700K  desktop: 1 to 29 in  38 seconds
i7-6820HQ laptop:  1 to 29 in  43 seconds
i5-3570K  desktop: 1 to 29 in 114 seconds

Puede hacer esto fácilmente de un solo subproceso y evitar esa dependencia de terceros utilizando el estándar para:

(for [start (range 0 end-idx step)]
  ... )

Bueno, el incorporado pmap incorporado también existe, pero Claypoole tiene más características y capacidad de ajuste.

NikoNyrh
fuente
Sí, hace que sea trivial distribuir. ¿Tendría tiempo de reevaluar mi solución? Estoy bastante seguro de que la obtendrá hasta 30 ahora. No tengo más optimizaciones a la vista.
NikoNyrh
Lamentablemente es un no por 30. Tiempo transcurrido: 217150.87386 ms
Ahaa, gracias por intentarlo: D Podría haber sido mejor ajustar una curva en esto e interpolar en qué valor decimal se gastan 120 segundos, pero aun así es un buen desafío.
NikoNyrh
1

Mathematica, n = 19

presione alt +. abortar y se imprimirá el resultado

k = 0;
For[n = 1, n < 1000, n++,
Z = Table[HammingDistance[#[[;; i]], #[[-i ;;]]], {i, Length@#}] & /@
Tuples[{0, 1}, n];
Table[If[Z[[i, j]] < 2, Z[[i, j]] = 0, Z[[i, j]] = 1], {i, 
Length@Z}, {j, n}];
k = Length@Union@Z]
Print["f(", n, ")=", k]
J42161217
fuente
No puedo ejecutar esto, ¿podría explicar cómo evita tomar tiempo exponencial? 2 ^ 241 es un número muy grande!
¿Puedes mostrar la salida del código?
1
Quise
1

Mathematica, 21

f [n_]: = Longitud @
     DeleteDuplicates @
      Transponer@
       Tabla [2> Tr @ IntegerDigits [#, 2] & / @ 
         BitXor [BitShiftRight [#, n - i], Mod [#, 2 ^ i]], {i, 1, 
         n - 1}] & @ Rango [0, 2 ^ (n - 1)];
Hacer [Imprimir [n -> f @ n], {n, Infinito}]

A modo de comparación, la respuesta de Jenny_mathy dan = 19 en mi equipo.

La parte más lenta es Tr@IntegerDigits[#, 2] &. Es una pena que Mathematica no tenga incorporado el peso de Hamming.


Si desea probar mi código, puede descargar una versión de prueba gratuita de Mathematica .

alephalpha
fuente
1

Versión de CA, utilizando una cuenta pop incorporada

Funciona mejor con clang -O3, pero también funciona si solo tienes gcc.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

unsigned long pairs(unsigned int n, unsigned long s) { 
  unsigned long result = 0;

  for(int d=1;d<=n;d++) { 
    unsigned long mx = 1 << d;
    unsigned long mask = mx - 1;

    unsigned long diff = (s >> (n - d)) ^ (s & mask);
    if (__builtin_popcountl(diff) <= 1)
      result |= mx;
  } 
  return result;

}

unsigned long f(unsigned long  n) { 
  unsigned long max = 1 << (n - 1);
#define BLEN (max / 2)
  unsigned char * buf = malloc(BLEN);
  memset(buf, 0, BLEN);
  unsigned long long * bufll = (void *) buf;

  for(unsigned long i=0;i<=max;i++) { 
    unsigned int r = pairs(n, i);
    buf[r / 8] |= 1 << (r % 8);
  } 

  unsigned long result = 0;

  for(unsigned long i=0;i<= max / 2 / sizeof(unsigned long long); i++) { 
    result += __builtin_popcountll(bufll[i]);
  } 

  free(buf);

  return result;
}

int main(int argc, char ** argv) { 
  unsigned int n = 1;

  while(1) { 
    printf("%d %ld\n", n, f(n));
    n++;
  } 
  return 0;
}
bartavelle
fuente
Llega a 24 muy rápidamente y luego termina. Necesitas aumentar el límite.
¡Dios mío, olvidé eliminar el código de referencia! Eliminaré las dos líneas ofensivas: /
bartavelle
@Lembik debería arreglarse ahora
bartavelle
1

Haskell, (no oficial n = 20)

Este es solo el enfoque ingenuo, hasta ahora sin optimizaciones. Me preguntaba qué tan bien le iría a otros idiomas.

Cómo usarlo (suponiendo que tenga instalada la plataforma haskell ):

  • Pegue el código en un archivo approx_corr.hs(o cualquier otro nombre, modifique los siguientes pasos en consecuencia)
  • Navega hasta el archivo y ejecuta ghc approx_corr.hs
  • correr approx_corr.exe
  • Ingrese el máximo n
  • Se muestra el resultado de cada cálculo, así como el tiempo real acumulado (en ms) hasta ese punto.

Código:

import Data.List
import Data.Time
import Data.Time.Clock.POSIX

num2bin :: Int -> Int -> [Int]
num2bin 0 _ = []
num2bin n k| k >= 2^(n-1) = 1 : num2bin (n-1)( k-2^(n-1))
           | otherwise  = 0: num2bin (n-1) k

genBinNum :: Int -> [[Int]]
genBinNum n = map (num2bin n) [0..2^n-1]

pairs :: [a] -> [([a],[a])]
pairs xs = zip (prefixes xs) (suffixes xs)
   where prefixes = tail . init . inits 
         suffixes = map reverse . prefixes . reverse 

hammingDist :: (Num b, Eq a) => ([a],[a]) -> b     
hammingDist (a,b) = sum $ zipWith (\u v -> if u /= v then 1 else 0) a b

f :: Int -> Int
f n = length $ nub $ map (map ((<=1).hammingDist) . pairs) $ genBinNum n
--f n = sum [1..n]

--time in milliseconds
getTime = getCurrentTime >>= pure . (1000*) . utcTimeToPOSIXSeconds >>= pure . round


main :: IO()
main = do 
    maxns <- getLine 
    let maxn = (read maxns)::Int
    t0 <- getTime 
    loop 1 maxn t0
     where loop n maxn t0|n==maxn = return ()
           loop n maxn t0
             = do 
                 putStrLn $ "fun eval: " ++ (show n) ++ ", " ++ (show $ (f n)) 
                 t <- getTime
                 putStrLn $ "time: " ++ show (t-t0); 
                 loop (n+1) maxn t0
falla
fuente
El código parece no dar salida mientras se ejecuta. Esto hace que sea un poco difícil de probar.
Extraño, ¿se compila sin error? ¿Qué sucede si intentas compilar el programa main = putStrLn "Hello World!"?
flawr
El Data.Bitsmódulo puede ser útil. Para su ciclo principal, podría usar algo como main = do maxn <- getmax; t0 <- gettime; loop 1dónde loop n|n==maxn = return ()y loop n = do printresult n (f n); t <- gettime; printtime (t-t0); loop (n+1). getmaxpodría, por ejemplo, usar getArgspara usar los argumentos del programa.
Christian Sievers
@ChristianSievers ¡Muchas gracias! Hice esta pregunta en stackoverflow, ¡creo que sería genial si pudieras agregar eso también!
flawr
No veo cómo responder allí. Ya tienes un bucle similar allí, y no dije nada sobre conseguir el tiempo: que ya tenías aquí.
Christian Sievers
1

Una solución de Haskell, usando popCount y paralelismo administrado manualmente

Compilar: ghc -rtsopts -threaded -O2 -fllvm -Wall foo.hs

(suelte el -llvmsi no funciona)

Correr : ./foo +RTS -N

module Main (main) where

import Data.Bits
import Data.Word
import Data.List
import qualified Data.IntSet as S 
import System.IO
import Control.Monad
import Control.Concurrent
import Control.Exception.Base (evaluate)

pairs' :: Int -> Word64 -> Int
pairs' n s = fromIntegral $ foldl' (.|.) 0 $ map mk [1..n]
  where mk d = let mask = 1 `shiftL` d - 1 
                   pc = popCount $! xor (s `shiftR` (n - d)) (s .&. mask)
               in  if pc <= 1 
                     then mask + 1 
                     else 0 

mkSet :: Int -> Word64 -> Word64 -> S.IntSet
mkSet n a b = S.fromList $ map (pairs' n) [a .. b]

f :: Int -> IO Int
f n 
   | n < 4 = return $ S.size $ mkSet n 0 mxbound
   | otherwise = do
        mvs <- replicateM 4 newEmptyMVar
        forM_ (zip mvs cpairs) $ \(mv,(mi,ma)) -> forkIO $ do
          evaluate (mkSet n mi ma) >>= putMVar mv
        set <- foldl' S.union S.empty <$> mapM readMVar mvs
        return $! S.size set
   where
     mxbound = 1 `shiftL` (n - 1)
     bounds = [0,1 `shiftL` (n - 3) .. mxbound]
     cpairs = zip bounds (drop 1 bounds)

main :: IO()
main = do
    hSetBuffering stdout LineBuffering
    mapM_ (f >=> print) [1..]
bartavelle
fuente
Parece que hay un problema de almacenamiento en búfer porque no obtengo ningún resultado si lo ejecuto desde la línea de comando cygwim.
Acabo de actualizar mi solución, pero no sé si ayudará mucho.
bartavelle
@Lembik No estoy seguro si eso es obvio, pero eso debería compilarse -O3y podría ser más rápido con -O3 -fllvm...
bartavelle
(Y todos los archivos de compilación deben eliminarse antes de volver a compilar, si no se
produce un
@Lembik: introduje el paralelismo. Debería ser un poco más rápido.
bartavelle
0

Python 2 + pypy, n = 22

Aquí hay una solución Python realmente simple como una especie de referencia de referencia.

import itertools
def hamming(A, B):
    n = len(A)
    assert(len(B) == n)
    return n-sum([A[i] == B[i] for i in xrange(n)])

def prefsufflist(P):
    n = len(P)
    return [hamming(P[:i], P[n-i:n]) for i in xrange(1,n+1)]

bound = 1
for n in xrange(1,25):
    booleans = set()
    for P in itertools.product([0,1], repeat = n):
        booleans.add(tuple(int(HD <= bound) for HD in prefsufflist(P)))
    print "n = ", n, len(booleans)

fuente