Проблемы с ленивой сверткой fn в Clojure

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

Рифтинг на версии jg-faustus, если бы я делал это в ветке эквивалента Clojure. Работает на меня. ~ 400 мс для 1 000 000 точек, ~ 25 мс для 100 000 на i7 Mackbook Pro.

Вероятная причина ошибок переполнения стека заключается в том, что ленивые преобразователи становятся слишком глубокими. (concat и map ленивы). Попробуйте обернуть эти вызовы в doall, чтобы принудительно оценить их возвращаемые значения.

Что касается более функционального решения, попробуйте что-то вроде этого:

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

Использование может использовать reduce для простого выполнения суммирования, а поскольку map создает ленивую последовательность, эта функция также ленива.

Не могу помочь с высокопроизводительной функциональной версией, но вы можете получить 100-кратное ускорение для императивной версии, отказавшись от лени и добавив подсказки типов:

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

С xs размером 100 КБ и is размером 2 ваш imp-convolve-1 занимает на моей машине ~ 6000 мс, когда он завернут в куклу, а этот - ~ 35 мс.

Обновить

Вот ленивая функциональная версия:

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

На размерах 100k и 2 он работает примерно 600 мс (в зависимости от 450-750 мс) против ~ 6000 мс для imp-convolve-1 и ~ 35 мс для imp-convolve-2.

Так что это функционально, лениво и имеет сносную производительность. Тем не менее, это в два раза больше кода, чем императивная версия, и мне потребовалось 1-2 дополнительных часа, чтобы найти, поэтому я не уверен, что вижу в этом смысл.

Я полностью за чистые функции, когда они делают код короче или проще или имеют какое-то другое преимущество перед императивной версией. Когда они этого не делают, я не возражаю переключиться в императивный режим.

Это одна из причин, почему я считаю, что Clojure великолепен, поскольку вы можете использовать любой подход по своему усмотрению.

Обновление 2:

Я внесу поправку в свой «какой смысл делать это функционально», сказав, что мне нравится эта функциональная реализация (вторая ниже на странице) Дэвида Кабаны.

Он короткий, удобочитаемый и занимает время ~140 мс при тех же размерах входных данных, что и выше (100 000 и 2), что делает его на сегодняшний день наиболее эффективной функциональной реализацией из тех, что я пробовал.

Учитывая, что он функциональный (но не ленивый), не использует подсказки типов и работает со всеми числовыми типами (а не только с двойными), это впечатляет.

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

Ленивое, функциональное решение можно выразить лаконично. Увы, производительность для > 2k нецелесообразна. Мне интересно посмотреть, есть ли способы ускорить его, не жертвуя читабельностью.

Изменить: после прочтения информативного сообщения drcabana на эту тему (http://erl.nfshost.com/2010/07/17/discrete-convolution-of-finite-vectors/), я обновил свой код, чтобы принимать векторы разных размеров. . Его реализация работает лучше: для размера xs 3 размер 1000000, ~ 2 с против ~ 3 с.

Правка 2: Взяв за основу идеи drcabana о простом изменении xs и padding is, я пришел к следующему:

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

Это, вероятно, настолько кратко, насколько это возможно, но в целом все еще медленнее, вероятно, из-за времени ожидания. Спасибо автору блога за взвешенный подход. Единственным преимуществом здесь является то, что приведенное выше действительно лениво в том смысле, что если я спрошу (nth res 10000), для получения результата потребуются только первые 10 000 вычислений.

На самом деле это не ответ ни на один из множества вопросов, которые вы задали, но у меня есть несколько комментариев к тем, которые вы не задавали.

1. Вероятно, вам не следует использовать nth для векторов. Да, это O(1), но поскольку nth работает с другими последовательностями за O(n), это (а) не дает понять, что вы ожидаете, что ввод будет вектором, и (б) означает, что если вы сделаете ошибка, ваш код таинственным образом станет очень медленным, а не сразу.

2. for и map оба ленивы, а aset - только побочные эффекты. Эта комбинация — верный путь к катастрофе: для побочных эффектов, подобных for, используйте doseq.

3. for и doseq допускают несколько привязок, поэтому вам не нужно накапливать их кучу, как вы (очевидно) делаете в Python.

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

   будет делать то, что вы хотите.

4. #(first %) более кратко записывается как first; аналогично #(+ %1 %2) равно +.

5. Вызов vec для набора промежуточных результатов, которым не обязательно быть векторами, замедлит работу. В частности, в vec-add достаточно вызывать vec только тогда, когда вы делаете окончательное возвращаемое значение: в (vec (reduce foo bar)) нет причин для foo превращать свои промежуточные результаты в векторы, если он никогда не использует их для произвольного доступа.