我正在写一些信号处理软件,我开始写一个
discrete convolution function.
这适用于前一万个左右的值列表,但随着它们变大(例如,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))))
解决方法
(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分支中这样做的话,重复j-g-faustus的版本.适合我.在i7 Mackbook Pro上,1,000,000点约为400毫秒,在100,000点上约为25毫秒.