在Clojure中使用延迟卷积fn的麻烦

Question

我正在编写一些信号处理软件,我开始编写一个离散卷积函数.

这适用于前一万个左右的值列表,但随着它们变大(例如,100k),我开始得到StackOverflow错误,当然.

不幸的是,我在将命令式卷积算法转换为递归和懒惰版本时遇到了很多麻烦,实际上使用速度足够快(至少有一点优雅也很好).

我也不是100%确定我完全没有这个功能,但是 - 如果我错过了什么/做错了什么,请告诉我.我认为这是正确的.

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

我将不得不寻求任何帮助:我仍然在Clojure中获得优势,并且使这个优雅,懒惰和/或递归会非常棒.

我有点惊讶的是,表达一种很容易用Lisp中的命令式语言表达的算法是多么困难.但也许我做错了!

编辑: 只是为了表明在命令式语言中表达是多么容易,并且为人们提供了运行良好且易于阅读的算法,这里是Python版本.除了更短,更简洁,更容易推理之外,它比Clojure代码执行的速度快几个数量级:甚至是使用Java数组的命令式Clojure代码.

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

另一方面,这是命令式Clojure代码.它还会从卷积中删除最后一个未完全浸入的值.因此除了缓慢和丑陋之外,它不是100%功能性的.也没有功能.

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

这太令人沮丧了.拜托,有人告诉我,我错过了一些明显的东西.

编辑3:

这是我昨天想到的另一个版本,展示了我希望如何能够表达它(虽然其他解决方案非常优雅;我只是在那里放另一个!)

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

它使用此实用程序功能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))))

Answer 1

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

如果我在 Clojure equiv 分支中执行此操作，则重复 jg-faustus 的版本。对我有用。在 i7 Mackbook Pro 上，1,000,000 个点约为 400 毫秒，100,000 个点约为 25 毫秒。