Functional programming – lazy convolution FN problem in clojure

I'm writing some signal processing software. I began to write a discrete revolution function

This applies to the first 10000 or so value lists, but as they get larger (for example, 100k), I start to get stackoverflow errors, of course

Unfortunately, I had a lot of trouble converting imperative convolution algorithms to recursion and amplifiers The lazy version is actually fast enough (at least a little elegant and good)

I'm not 100% sure I don't have this feature at all, but – if I miss / do something wrong, please let me know I think this is correct

(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)))))))

I will have to ask for any help: I still gain an advantage in clojure, and it would be great to make this elegant, lazy and / or recursive

I'm a little surprised how difficult it is to express an algorithm that is easy to express in the imperative language in lisp But maybe I did it wrong!

Editor: just to show how easy it is to express in imperative language and provide people with a well run and easy to read algorithm. Here is the python version In addition to being shorter, simpler, and easier to reason, it executes several orders of magnitude faster than clojure Code: even imperative clojure code using java arrays

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

On the other hand, this is imperative clojure code It also removes the last incompletely immersed value from the convolution So it's not 100% functional except slow and ugly It has no function

(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)))))))))

This is so frustrating Please, someone told me I missed something obvious

Edit 3:

This is another version I thought of yesterday, showing how I hope to express it (although other solutions are very elegant; I just put another one there!)

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

It uses this utility 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))))

Solution

(defn ^{:static true} convolve ^doubles [^doubles xs ^doubles is]

(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))

If I do this in the clojure equiv branch, repeat the version of j-g-faustus Fit me. On i7 mackbook pro, 400 ms at 1000000 points and 25 ms at 100000 points