;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; $Id: mnist-kernel.lsh,v 1.9 2005/05/19 15:30:09 leonb Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(libload "svm/kernel")
(libload "libidx/idx-float")
(libload "libidx/idx-double")
(libload "libidx/idx-int")

(defclass MnistRBFKernel SVMKernel
  ((-idx2- (-ubyte-)) x)
  ((-idx1- (-ubyte-)) y)
  ((-double-) gamma) 
  ((-int-) count) )

(defmethod MnistRBFKernel MnistRBFKernel(trx try tex tey g)
  ((-idx3- (-ubyte-)) trx)
  ((-idx1- (-ubyte-)) try)
  ((-idx3- (-ubyte-)) tex)
  ((-idx1- (-ubyte-)) tey)
  ((-double-) g) 
  (setq x (ubyte-matrix 70000 (* 28 28)))
  (setq y (ubyte-matrix 70000))
  (copy-any-matrix trx (narrow x 0 60000 0))
  (copy-any-matrix tex (narrow x 0 10000 60000))
  (copy-any-matrix try (narrow y 0 60000 0))
  (copy-any-matrix tey (narrow y 0 10000 60000))
  (setq gamma g)
  (setq count 0) )

(defmethod MnistRBFKernel vector(i)
  ((-int-) i)
  (select x 0 i) )

(defmethod MnistRBFKernel label(i)
  ((-int-) i)
  (to-double (y i)) )

(defmethod MnistRBFKernel call(i j)
  ((-int-) i j)
  (incr count)
  (let ((sum 0)
	(x1 (select x 0 i))
	(x2 (select x 0 j)) )
    ((-int-) sum)
    #{{ unsigned char *p1 =  IDX_PTR ($x1, unsigned char);
        unsigned char *p2 =  IDX_PTR ($x2, unsigned char);
        $sum = asm_sumsqr_784(p1,p2); }#}
    (qexpmx (* (/ gamma 16384) sum)) ) )

(delete mnist-kernel)

(let ((dhc-make-lushflags (concat dhc-make-lushflags " -g -funroll-all-loops")))
  (dhc-make () 
            #{
                static int 
                asm_sumsqr_784(unsigned char *p1, unsigned char *p2)
                {
                  int i;
                  int s=0;

                #if defined(__SSE2__)

                  asm("pxor %%xmm6,%%xmm6\n\t"
                      "pxor %%xmm7,%%xmm7"
                      : : : "memory" );
                  for(i=0; i<784; i+=16)
                    {
                      asm("movdqa %0,%%xmm0\n\t"
                          "movdqa %1,%%xmm1\n\t"
                          "prefetcht1 %2\n\t"
                          "prefetcht1 %3\n\t" 
                          "movdqa %%xmm0,%%xmm2\n\t"
                          "movdqa %%xmm1,%%xmm3\n\t"
                          "punpcklbw %%xmm7,%%xmm2\n\t" 
                          "punpcklbw %%xmm7,%%xmm3\n\t" 
                          "punpckhbw %%xmm7,%%xmm0\n\t" 
                          "punpckhbw %%xmm7,%%xmm1\n\t"
                          "psubw %%xmm2,%%xmm3\n\t"
                          "psubw %%xmm0,%%xmm1\n\t"
                          "pmaddwd %%xmm3,%%xmm3\n\t"
                          "pmaddwd %%xmm1,%%xmm1\n\t"
                          "paddd %%xmm3,%%xmm6\n\t"
                          "paddd %%xmm1,%%xmm6"
                          : : "m"(p1[i]), "m"(p2[i]), 
                              "m"(p1[i+64]), "m"(p2[i+64]) );
                    }
                  asm("pshufd %1,%%xmm6,%%xmm7\n\t"
                      "paddd %%xmm7,%%xmm6\n\t"
                      "pshufd %2,%%xmm6,%%xmm7\n\t"
                      "paddd %%xmm7,%%xmm6\n\t"
                      "movd  %%xmm6,%0"
                      : "=m" (s) : "i"(0x4e), "i"(0x93) : "memory" );

                #elif defined(__MMX__)

                  asm("pxor %%mm6,%%mm6\n\t"
                      "pxor %%mm7,%%mm7"
                      : : : "memory" );
                  for(i=0; i<784; i+=8)
                    {
                      asm("movq %0,%%mm0\n\t"
                          "movq %1,%%mm1\n\t"
                # if defined(__SSE__)
                          "prefetcht1 %2\n\t"
                          "prefetcht1 %3\n\t" 
                # endif
                          "movq %%mm0,%%mm2\n\t"
                          "movq %%mm1,%%mm3\n\t"
                          "punpcklbw %%mm7,%%mm2\n\t" 
                          "punpcklbw %%mm7,%%mm3\n\t" 
                          "punpckhbw %%mm7,%%mm0\n\t" 
                          "punpckhbw %%mm7,%%mm1\n\t"
                          "psubw %%mm2,%%mm3\n\t"
                          "psubw %%mm0,%%mm1\n\t"
                          "pmaddwd %%mm3,%%mm3\n\t"
                          "pmaddwd %%mm1,%%mm1\n\t"
                          "paddd %%mm3,%%mm6\n\t"
                          "paddd %%mm1,%%mm6"
                          : : "m"(p1[i]), "m"(p2[i])
                # if defined(__SSE__)
                            , "m"(p1[i+64]), "m"(p2[i+64]) 
                # endif
                          );
                    }
                  asm("movq  %%mm6,%%mm7\n\t"
                      "psrlq %1,%%mm7\n\t"
                      "paddd %%mm6,%%mm7\n\t"
                      "movd  %%mm7,%0\n\t"
                      "emms"
                      : "=m" (s) : "i"(32) : "memory" );
                #else

                  for(i=0; i<784; i++)
                    {
                      int d = (int)p1[i] - (int)p2[i];
                      s += d*d;
                    }

                #endif
                  return s;
                }
  #}
  (MnistRBFKernel MnistRBFKernel vector label call) ) )

(defvar mnist-train-x (load-matrix "mnist/train-images-idx3-ubyte"))
(defvar mnist-test-x (load-matrix "mnist/t10k-images-idx3-ubyte")) 
(defvar mnist-train-y (load-matrix "mnist/train-labels-idx1-ubyte")) 
(defvar mnist-test-y (load-matrix "mnist/t10k-labels-idx1-ubyte"))
(defvar mnist-kernel ())

(de load-mnist() 
  (delete mnist-kernel)
  (setq mnist-kernel 
        (new MnistRBFKernel 
             mnist-train-x mnist-train-y
             mnist-test-x mnist-test-y
             0.005 ) ) )

(load-mnist)
