6

Estoy escribiendo un software de procesamiento de señal, y estoy empezando por escribir un discrete convolution function.Problemas con lazy convolution fn en Clojure

Esto funciona bien para las primeras diez mil listas de valores, pero a medida que crecen (digamos, 100k), empiezo a tener errores de StackOverflow, por supuesto.

Por desgracia, estoy teniendo un montón de problemas para convertir el algoritmo de convolución imperativo que tengo a un recursiva & versión perezosa que en realidad es lo suficientemente rápido para usar (que tiene al menos un mínimo de elegancia sería bueno también)

Tampoco estoy 100% seguro de tener esta función completamente correcta, pero – por favor avíseme si me falta algo/estoy haciendo algo mal. I piensa es correcto.

(defn convolve 
    " 
    Convolves xs with is. 

    This is a discrete convolution. 

    'xs :: list of numbers 
    'is :: list of numbers 
    " 
    [xs is] 
    (loop [xs xs finalacc() acc()] 
    (if (empty? xs) 
     (concat finalacc acc) 
     (recur (rest xs) 
      (if (empty? acc) 
       () 
       (concat finalacc [(first acc)])) 
      (if (empty? acc) 
       (map #(* (first xs) %) is) 
       (vec-add 
       (map #(* (first xs) %) is) 
       (rest acc))))))) 

Yo estaría muy agradecido por cualquier tipo de ayuda: Todavía estoy orientarme en Clojure, y haciendo de este elegante y perezoso y/o recursiva sería maravilloso.

Estoy un poco sorprendido de lo difícil que es expresar un algoritmo que es bastante fácil de expresar en un lenguaje imperativo en un Lisp. ¡Pero quizás lo estoy haciendo mal!

EDIT: sólo para mostrar lo fácil que es para expresar en un lenguaje imperativo, y para dar a la gente el algoritmo que funciona muy bien y es fácil de leer, aquí está la versión de Python. Además de ser más corto, más conciso y mucho más fácil de razonar, ejecuta órdenes de magnitud más rápido que el código Clojure: incluso mi código imperativo Clojure con matrices Java.

from itertools import repeat 

def convolve(ns, ms): 
    y = [i for i in repeat(0, len(ns)+len(ms)-1)] 
    for n in range(len(ns)): 
     for m in range(len(ms)): 
      y[n+m] = y[n+m] + ns[n]*ms[m] 
    return y 

Aquí, por otro lado, está el código imperativo Clojure. También elimina los últimos valores, no completamente sumergidos, de la convolución. Además de ser lento y feo, no es 100% funcional. Ni funcional

(defn imp-convolve-1 
    [xs is] 
    (let [ys (into-array Double (repeat (dec (+ (count xs) (count is))) 0.0)) 
     xs (vec xs) 
     is (vec is)] 
    (map #(first %) 
      (for [i (range (count xs))] 
      (for [j (range (count is))] 
       (aset ys (+ i j) 
        (+ (* (nth xs i) (nth is j)) 
         (nth ys (+ i j))))))))) 

Esto es tan desalentador. Por favor, que alguien me enseñe que me he perdido algo obvio.

Datos 3:

Aquí hay otra versión que pensé ayer, mostrando cómo me gustaría ser capaz de expresarlo (aunque otras soluciones son muy elegantes, yo sólo voy a poner otra por ahí!)

(defn convolve-2 
    [xs is] 
    (reduce #(vec-add %1 (pad-l %2 (inc (count %1)))) 
     (for [x xs] 
      (for [i is] 
      (* x i))))) 

Se utiliza esta función de utilidad vec-add:

(defn vec-add 
    ([xs] (vec-add xs [])) 
    ([xs ys] 
    (let [lxs (count xs) 
      lys (count ys) 
      xs (pad-r xs lys) 
      ys (pad-r ys lxs)] 
     (vec (map #(+ %1 %2) xs ys)))) 
    ([xs ys & more] 
    (vec (reduce vec-add (vec-add xs ys) more)))) 
    (vec (reduce vec-add (vec-add xs ys) more)))) 
+0

Esto no es en absoluto una respuesta a su pregunta, pero ... la implementación de la convolución mediante el uso de transformadas rápidas de Fourier reduce considerablemente el número de operaciones (consulte http://stackoverflow.com/questions/3084987/low- pass-filter-using-fft-instead-of-convolution-implementation, por ejemplo). – mtrw

+0

No sé lo suficiente sobre los lenguajes funcionales para saber si esto ayudará, pero podría estar interesado en http://stackoverflow.com/questions/803055/how-do-i-do-convolution-in-f. – mtrw

+0

@mtrw En producción, esto definitivamente se hará usando la FFT; por ahora, sin embargo, me encantaría tener una solución agradable y funcional. Pero no estoy seguro de lo posible que es (en Clojure). – Isaac

Respuesta

4
(defn ^{:static true} convolve ^doubles [^doubles xs ^doubles is] 
    (let [xlen (count xs) 
     ilen (count is) 
     ys (double-array (dec (+ xlen ilen)))] 
    (dotimes [p xlen] 
     (dotimes [q ilen] 
     (let [n (+ p q), x (aget xs p), i (aget is q), y (aget ys n)] 
      (aset ys n (+ (* x i) y))))) 
    ys)) 

Riffing en la versión de j-g-faustus si estaba haciendo esto en la rama equiv Clojure. Funciona para mi. ~ 400ms para 1,000,000 puntos, ~ 25ms para 100,000 en un i7 Mackbook Pro.

4

La causa probable de la XX Los errores de desbordamiento de la pila es que los thigs perezosos se están volviendo demasiado profundos. (concat y map son flojos). Intente envolver esas llamadas en doall para forzar la evaluación de sus valores devueltos.

En cuanto a una solución más funcional, intentar algo como esto:

(defn circular-convolve 
    "Perform a circular convolution of vectors f and g" 
    [f g] 
    (letfn [(point-mul [m n] 
     (* (f m) (g (mod (- n m) (count g))))) 
     (value-at [n] 
     (reduce + (map #(point-mul % n) (range (count g)))))] 
    (map value-at (range (count g))))) 

Uso puede utilizar para llevar a cabo reduce suma facilidad, y desde map produce una secuencia perezosa, esta función también es perezoso.

+0

Desafortunadamente, envolver mis 'maps' y' concats' en 'doalls' hace que transcurra un tiempo inaceptablemente prolongado (¡no estoy seguro de que termine!) Para calcular una convolución con una matriz con 100.000 puntos y la otra con 2. .. Su circunvolución circular no devuelve listas de la longitud correcta o valores correctos (¡no lo creo!) ¡Gracias por la ayuda! – Isaac

4

no puede ayudar con una versión funcional de alto rendimiento, pero se puede conseguir un aumento de velocidad de 100 veces para la versión imperativa por lo que dejaría la pereza y añadiendo toques de tipo:

(defn imp-convolve-2 [xs is] 
    (let [^doubles xs (into-array Double/TYPE xs) 
     ^doubles is (into-array Double/TYPE is) 
     ys (double-array (dec (+ (count xs) (count is)))) ] 
    (dotimes [i (alength xs)] 
     (dotimes [j (alength is)] 
     (aset ys (+ i j) 
      (+ (* (aget xs i) (aget is j)) 
      (aget ys (+ i j)))))) 
    ys)) 

Con xs 100k tamaño y is tamaño 2, su imp-convolve-1 toma ~ 6,000ms en mi máquina cuando está envuelta en una doall, mientras que esta toma ~ 35ms.

Actualizar

He aquí una versión funcional perezoso:

(defn convolve 
    ([xs is] (convolve xs is [])) 
    ([xs is parts] 
    (cond 
     (and (empty? xs) (empty? parts)) nil 
     (empty? xs) (cons 
        (reduce + (map first parts)) 
        (convolve xs is 
         (remove empty? (map rest parts)))) 
     :else (cons 
       (+ (* (first xs) (first is)) 
       (reduce + (map first parts))) 
       (lazy-seq 
       (convolve (rest xs) is 
        (cons 
        (map (partial * (first xs)) (rest is)) 
        (remove empty? (map rest parts))))))))) 

En los tamaños 100 mil y 2, que los relojes en ~ 600 ms (variando 450-750ms) vs ~ 6,000ms para Imp convolve-1 y ~ 35ms para imp-convolve-2.

Por lo tanto, es funcional, vago y tiene un rendimiento tolerable. Aún así, es el doble de código que la versión imperativa y me tomó 1-2 horas adicionales para encontrar, así que no estoy seguro de que vea el punto.

Soy todo para funciones puras cuando hacen el código más corto o más simple, o tienen algún otro beneficio sobre una versión imperativa. Cuando no lo hacen, no tengo inconveniente en cambiar al modo imperativo.

Cuál es una de las razones por las que creo que Clojure es genial, ya que puede usar cualquiera de los enfoques que considere oportuno.

Actualización 2:

voy a enmendar mi "cuál es el punto de hacer esto funcionalmente" diciendo que me gusta this functional implementation (la segunda, final de la página) de David Cabana.

Es breve, legible y tiene un tiempo de ~ 140ms con los mismos tamaños de entrada que el anterior (100,000 y 2), por lo que es la implementación funcional con mejor rendimiento de los que he probado.

Teniendo en cuenta que es funcional (pero no flojo), no utiliza sugerencias de tipo y funciona para todos los tipos numéricos (no solo dobles), eso es bastante impresionante.

+0

¡Lo aprecio, prácticamente no tengo experiencia haciendo Clojure imperativo! Tal vez debería aprender ... ¡Gracias! (No estoy aceptando, ya que sigo buscando una versión funcional, pero tengo mi +1) – Isaac

+0

Aunque, curiosamente, cuando lo intento con 1 millón de puntos, aparece un "error en el filtro del proceso: Número incorrecto de argumentos: nada, 0" . Aunque no hay errores de Java que pueda ver ... Funciona bien para 100k, sin embargo. La versión de My Python funciona para 1 millón; no la he probado con más. – Isaac

+0

Debería funcionar para puntos enteros/MAX_VALUE (aproximadamente 2 mil millones), es el tamaño máximo de una matriz de Java. No veo nada diferente por 1 millón de puntos. Claro que no sucede en otro lado? "Error en el filtro de proceso" no suena como si proviniera de la función convolve. –

3
(defn convolve [xs is] 
    (if (> (count xs) (count is)) 
    (convolve is xs) 
    (let [os (dec (+ (count xs) (count is))) 
      lxs (count xs) 
      lis (count is)] 
     (for [i (range os)] 
     (let [[start-x end-x] [(- lxs (min lxs (- os i))) (min i (dec lxs))] 
       [start-i end-i] [(- lis (min lis (- os i))) (min i (dec lis))]] 
      (reduce + (map * 
         (rseq (subvec xs start-x (inc end-x))) 
         (subvec is start-i (inc end-i))))))))) 

Es posible expresar una solución perezosa y funcional en términos concisos. Por desgracia, el rendimiento de> 2k no es práctico. Me interesa ver si hay formas de acelerarlo sin sacrificar la legibilidad.

Editar: Después de leer el post informativo de drcabana sobre el tema (http://erl.nfshost.com/2010/07/17/discrete-convolution-of-finite-vectors/), He actualizado mi código para aceptar diferentes vectores de tamaño.Su aplicación es más funcional: para el tamaño XS 3, es el tamaño de 1000000, ~ 2s vs ~ 3s

Edición 2: Tomando las ideas de drcabana de simplemente invirtiendo XS y el relleno es, yo llegamos a:

(defn convolve [xs is] 
    (if (> (count xs) (count is)) 
    (convolve is xs) 
    (let [is (concat (repeat (dec (count xs)) 0) is)] 
     (for [s (take-while not-empty (iterate rest is))] 
     (reduce + (map * (rseq xs) s)))))) 

Esto es probablemente tan conciso como se va a obtener, pero aún es más lento en general, probablemente debido a la toma de tiempo. Felicitaciones al autor del blog por un enfoque bien considerado. La única ventaja aquí es que lo anterior es realmente flojo porque si pregunto (enésimo res 10000), solo necesitará los primeros cálculos de 10k para llegar a un resultado.

3

No es realmente una respuesta a ninguna de las muchas preguntas que ha formulado, pero tengo varios comentarios sobre las que no ha preguntado.

  1. Probablemente no deba usar nth en vectores. Sí, es O (1), pero como nth funciona en otras secuencias en O (n), (a) no deja en claro que espera que la entrada sea un vector, y (b) significa que si hace una error, su código misteriosamente se volverá realmente lento en lugar de fallar inmediatamente.

  2. for y map son flojos, y aset son efectos secundarios solamente. Esta combinación es una receta para el desastre: para comportamientos similares al for, utilice doseq.

  3. for y doseq permiten múltiples enlaces, por lo que no necesita acumular montones de ellos como usted (al parecer) hacer en Python.

    (doseq [i (range (count cs)) 
         j (range (count is))] 
        ...) 
    

    hará lo que quiera.

  4. #(first %) es más concisamente escrito como first; del mismo modo #(+ %1 %2) es +.

  5. Calling vec en un montón de resultados intermedios que no necesidad de ser vectores le ralentizarán. Específicamente en vec-add es suficiente llamar solo al vec cuando se realiza un valor de retorno final: en (vec (reduce foo bar)) no hay razón para que foo convierta sus resultados intermedios en vectores si nunca los usa para acceso aleatorio.

Cuestiones relacionadas