clyce/apcluster-whatswrong.clj

## apcluster-whatswrong.clj
(require '[clojure.core.matrix :refer :all])

(set-current-implementation :vectorz)

(defn euclidien-sqr [x y]
  (-
    (apply +
           (map #(* % %)
                (sub! (array x) (array y))))))

(defn m-median [m]
  (let [vs (sort (apply concat m))
        n (inc (count vs))]
    (if (pos? (rem n 2))
      (let [t (quot n 2)]
        (/ (+ (nth vs (dec t)) (nth vs t)) 2))
      (nth vs (dec (/ n 2))))))

(defn v-max
  [v]
  (loop [n 0
         x (first v)
         v (rest v)
         res-idx n
         res-x x]
    (let [bigger? (> x res-x)
          res-idx (if bigger? n res-idx)
          res-x (if bigger? x res-x)]
      (if (empty? v)
        [res-idx res-x]
        (recur
          (inc n) (first v) (rest v)
          res-idx res-x)))))

(defn m-max
  [m d]
  (mapv v-max
        (case d 1 (rows m) 0 (columns m))))

(defn ap-clustering*
  [coll
   & {:keys [dist-fn convergence-iter max-iter damping]
      :or {dist-fn euclidien-sqr convergence-iter 15 max-iter 200 damping 0.5}}]
  (let [n (count coll)
        s (new-array [n n])
        idxs (range n)
        _ (doseq [i idxs]
            (let [x (nth coll i)]
              (set-row!
                s i
                (array
                  (map (partial dist-fn x) coll)))))
        preference (m-median s)
        _ (doseq [x idxs]
            (mset! s x x preference))]
    (println "s" s)
    (loop [a (new-array [n n])
           r (new-array [n n])
           centers #{}
           current-convergence-iter 0
           iter 0]
      (let [;;calc r
            tmp (add a s)
            _ (println "init: " tmp)

            [I Y] ((juxt (partial map first)
                         (partial map second))
                   (m-max tmp 1))

            _ (doseq [i idxs]
                (mset! tmp i (nth I i)
                       (- Double/MAX_VALUE)))

            [I2 Y2] ((juxt (partial map first)
                           (partial map second))
                     (m-max tmp 1))


            tmp (sub s (transpose (array (repeat n Y))))

            _ (doseq [i idxs]
                (let [j (nth I i)]
                  (mset! tmp i j
                         (- (mget s i j) (nth Y2 i)))))

            _ (println "r" tmp)

            r (add (mul damping r)
                   (mul (- 1 damping) tmp))
            _ (println "real-r" r)

            ;;calc a
            tmp (new-array [n n])

            _ (doseq [i idxs j idxs]
                (mset! tmp i j
                       (if (= i j)
                         (mget r i j)
                         (max 0
                              (mget r i j)))))

            _ (println "a*" tmp)

            _ (println "tmp"(array (repeat n (map (partial apply +) (columns tmp)))))
            tmp (sub tmp (array (repeat n (map (partial apply +) (columns tmp)))))
            _ (println "a" tmp)

            da (diagonal tmp)

            _ (doseq [i idxs]
                (mset! tmp i i (mget da i)))

            _ (println "aa" tmp)

            a (sub (mul damping a)
                   (mul (- 1 damping) tmp))

            _ (println "real-a" a)
            t (diagonal (add r a))

            _ (println "t" t)

            new-centers (set (filter #(pos? (mget t %)) idxs))]

        (if (and (< iter max-iter)
                 (< current-convergence-iter convergence-iter))
          (recur a r
                 new-centers
                 (if (= centers new-centers)
                   (inc current-convergence-iter)
                   0)
                 (inc iter))
          (for [i idxs]
            [(nth coll i)
             (if (centers i)
               i
               (first
                 (last
                   (sort-by last
                            (map (fn [k] [k (mget a i k)]) centers)))))]))))))
	(require '[clojure.core.matrix :refer :all])

	(set-current-implementation :vectorz)

	(defn euclidien-sqr [x y]
	(-
	(apply +
	(map #(* % %)
	(sub! (array x) (array y))))))

	(defn m-median [m]
	(let [vs (sort (apply concat m))
	n (inc (count vs))]
	(if (pos? (rem n 2))
	(let [t (quot n 2)]
	(/ (+ (nth vs (dec t)) (nth vs t)) 2))
	(nth vs (dec (/ n 2))))))

	(defn v-max
	[v]
	(loop [n 0
	x (first v)
	v (rest v)
	res-idx n
	res-x x]
	(let [bigger? (> x res-x)
	res-idx (if bigger? n res-idx)
	res-x (if bigger? x res-x)]
	(if (empty? v)
	[res-idx res-x]
	(recur
	(inc n) (first v) (rest v)
	res-idx res-x)))))

	(defn m-max
	[m d]
	(mapv v-max
	(case d 1 (rows m) 0 (columns m))))

	(defn ap-clustering*
	[coll
	& {:keys [dist-fn convergence-iter max-iter damping]
	:or {dist-fn euclidien-sqr convergence-iter 15 max-iter 200 damping 0.5}}]
	(let [n (count coll)
	s (new-array [n n])
	idxs (range n)
	_ (doseq [i idxs]
	(let [x (nth coll i)]
	(set-row!
	s i
	(array
	(map (partial dist-fn x) coll)))))
	preference (m-median s)
	_ (doseq [x idxs]
	(mset! s x x preference))]
	(println "s" s)
	(loop [a (new-array [n n])
	r (new-array [n n])
	centers #{}
	current-convergence-iter 0
	iter 0]
	(let [;;calc r
	tmp (add a s)
	_ (println "init: " tmp)

	[I Y] ((juxt (partial map first)
	(partial map second))
	(m-max tmp 1))

	_ (doseq [i idxs]
	(mset! tmp i (nth I i)
	(- Double/MAX_VALUE)))

	[I2 Y2] ((juxt (partial map first)
	(partial map second))
	(m-max tmp 1))


	tmp (sub s (transpose (array (repeat n Y))))

	_ (doseq [i idxs]
	(let [j (nth I i)]
	(mset! tmp i j
	(- (mget s i j) (nth Y2 i)))))

	_ (println "r" tmp)

	r (add (mul damping r)
	(mul (- 1 damping) tmp))
	_ (println "real-r" r)

	;;calc a
	tmp (new-array [n n])

	_ (doseq [i idxs j idxs]
	(mset! tmp i j
	(if (= i j)
	(mget r i j)
	(max 0
	(mget r i j)))))

	_ (println "a*" tmp)

	_ (println "tmp"(array (repeat n (map (partial apply +) (columns tmp)))))
	tmp (sub tmp (array (repeat n (map (partial apply +) (columns tmp)))))
	_ (println "a" tmp)

	da (diagonal tmp)

	_ (doseq [i idxs]
	(mset! tmp i i (mget da i)))

	_ (println "aa" tmp)

	a (sub (mul damping a)
	(mul (- 1 damping) tmp))

	_ (println "real-a" a)
	t (diagonal (add r a))

	_ (println "t" t)

	new-centers (set (filter #(pos? (mget t %)) idxs))]

	(if (and (< iter max-iter)
	(< current-convergence-iter convergence-iter))
	(recur a r
	new-centers
	(if (= centers new-centers)
	(inc current-convergence-iter)
	0)
	(inc iter))
	(for [i idxs]
	[(nth coll i)
	(if (centers i)
	i
	(first
	(last
	(sort-by last
	(map (fn [k] [k (mget a i k)]) centers)))))]))))))