incanter.contrib.vecmath.mvector

mvector.clj

(ns
  ^{:doc
    "Implementation of vector products using DoubleMatrix1D from Parallel Colt"}
  incanter.contrib.vecmath.mvector
  (:use [incanter.core :only (matrix
                              ncol nrow dim
                              plus minus mult div
                              mmult trans)])
  (:import [cern.colt.matrix.tdouble.algo SmpDoubleBlas]
           [cern.colt.matrix.tdouble DoubleMatrix1D DoubleMatrix2D]
           [cern.colt.matrix.tdouble.impl DenseDoubleMatrix1D]
           [incanter Matrix]))



(def ^:dynamic *blas* (SmpDoubleBlas.))

(defn mvector
  "Returns an instance of 1D Colt vector (DoubleMatrix1D)."
  [xs]
  (cond
    (instance? DoubleMatrix1D xs)
    xs ; unchanged
    (and (instance? DoubleMatrix2D xs)
         (some #(= 1 %) (dim xs)))
    (.vectorize xs) ;
    (sequential? xs)
    (DenseDoubleMatrix1D. (into-array Double/TYPE xs))
    :else
    (throw (IllegalArgumentException.
            "not a single-row or single-column matrix nor a sequence"))))

(defn as-matrix
  "Convert to Incanter matrix if necessary."
  [x]
  (cond
    (instance? Matrix x) x
    (instance? DoubleMatrix1D x) (matrix (.toArray x))
    :else (matrix x)))

(defn dot
  "Dot product between two vectors.

  Vectors can be 1D Colt vectors (DoubleMatrix1D), single-row or
  single-column 2D matrices (DoubleMatrix2D) or sequences.
  "
  [u v]
  (letfn [(mmult-dot [u v]
            "Dot product between single-column 2D matrices
             (using Incanter's mmult)."
            (mmult (trans u) v))
          (to-column [mat]
            (cond
              (= 1 (ncol mat)) mat
              (= 1 (nrow mat)) (trans mat)
              :else (throw (IllegalArgumentException.
                            "matrix is not single row or single-column"))))]
    (cond
      (every? #(instance? DoubleMatrix1D %) [u v])
      (.ddot *blas* u v)
      (every? #(instance? DoubleMatrix2D %) [u v])
      (mmult-dot (to-column u) (to-column v))
      :else
      (let [u' (mvector u)
            v' (mvector v)]
        (.ddot *blas* u' v')))))

(defn cross
  "Cross product between two vectors.

  Implementation:
  http://en.wikipedia.org/wiki/Cross_product#Conversion_to_matrix_multiplication
  "
  [u v]
  (letfn [(to-seq [u] (if (instance? DoubleMatrix1D u) (.toArray u) u))]
    (let [u' (to-seq u)
          v' (to-seq v)
          [a1 a2 a3] u'
          m (matrix [[0 (- a3) a2]
                     [a3 0 (- a1)]
                     [(- a2) a1 0]])]
      (mmult m (matrix v')))))