;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; 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: idx-convol.lsh,v 1.3 2002/10/09 06:43:39 profshadoko Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? **** Convolution, Subsampling, Oversampling.
;; This library contains functions for performing 1D and 2D convolutions,
;; backconvolutions, subsampling with local averaging, and 
;; oversampling by sample replication.
;; Each functionality is provided as a type-independent
;; splicing macro (DMD) with names of the form <midx-m1xxx>
;; and <midx-m2xxx>, as well as type specific functions
;; for floats, doubles, and bytes (with names of the form
;; <idx-f1xxx>, <idx-d1xxx>, <idx-u1xxx>.
;; Two kinds of convolution functions are provided: generic
;; 1D and 2D convolutions for any size kernel, and optimized
;; 1D convolutions with kernels of size 2,3,4 and 5.
;; The optimized versions are faster than the regular ones
;; by about a factor of 2 through the use of loop unrolling 
;; and circular register stacks.

#? *** Fast 1D Convolutions
;; a library of 1D convolutions that uses automatically generated
;; code with unrolled loops and circular register stacks
;; to limit memory accesses and speed up execution.



#? (midx-m1fastconvolacc <in> <kernel> <out> <ksize> <type>)
;; ultrafast 1D convolution with accumulation of result into the output.
;; <in> must be a CONTIGUOUS idx1, <kernel> an idx1, <out> a CONTIGUOUS idx1,
;; <ksize> must be the size of the kernel, and <type> a string
;; containing a valid C type, e.g. "float" (preferably the type
;; of the objects in the 3 vectors).
;; This macro is expanded at compile time into a bunch of "cinline"
;; instructions. It unrolls loops with a circular register stack
;; to minimize index calculations and memory access.
;; Though this code is particularly fast on processors that
;; have lots of floating point registers (MIPS, etc), it manages to
;; to beat the pedestrian version by about a factor of 2 on pentiums
;; (reaching approx 200MFLOPS per GHz for kernel of size 5).
;; THIS CODE DOES NOT DO ANY RUN-TIME CHECK FOR ARRAY SIZE, SO BE CAREFULL.

(dmd midx-m1fastconvolacc (in kernel out ksize type)
  (cons `progn 
	(append
	 `((cinline ,(sprintf "{ %s *r, *p, *ker;" type)))
	 (all ((x (range 0 (1- ksize))))
	   `(cinline ,(sprintf "  register %s w%d, x%d;" type x x)))
	 `((cinline " register int i;")
	   (cinline " int kmod, outdim, n, nmodk;")
	   (cinline "kmod = (%s)->mod[0];" ,kernel)
	   (cinline "outdim = (%s)->dim[0];" ,out)
	   (cinline ,(sprintf "n = outdim / %d;" ksize))
	   (cinline ,(sprintf "nmodk = outdim %%%% %d;" ksize))
	   (cinline ,(concat "  ker = IDX_PTR((%s)," type ");") ,kernel)
	   (cinline ,(concat "  p   = IDX_PTR((%s)," type ");") ,in)
	   (cinline ,(concat "  r   = IDX_PTR((%s)," type ");") ,out))
	 ;; initialize weights in registers
	 (all ((x (range 0 (1- ksize))))
	   `(cinline ,(sprintf "  w%d = *ker; ker += kmod;" x)))
	 ;; initialize inputs in registers
	 (all ((x (range 0 (- ksize 2))))
	   `(cinline ,(sprintf "  x%d = *p++;" x)))
	 `((cinline "for (i = 0; i<n; i++) {"))
	 (let ((lx (range 0 (1- ksize)))
	       (lw (range 0 (1- ksize)))
	       (r ()))
	   (repeat ksize
	     (setq r (nconc1 r `(cinline ,(sprintf "  x%d = *p++;" (last lx)))))
	     (setq r (nconc1 
		      r 
		      `(cinline 
			,(concat 
			  (sprintf "*r++ += w%d * x%d " (car lw) (car lx))
			  (apply concat 
				 (all ((w (cdr lw)) (x (cdr lx)))
				   (sprintf " + w%d * x%d " w x)))
			  ";"))))
	     (setq lx (cdr (nconc1 lx (car lx)))))
	   r)
	 `((cinline "}"))
	 ;; add the remainder of the convolutions
	 (let ((lx (range 0 (1- ksize)))
	       (lw (range 0 (1- ksize)))
	       (r ()))
	   (repeat (1- ksize)
	     (setq r (nconc1 r `(cinline "if (nmodk-- == 0) goto end;")))
	     (setq r (nconc1 r `(cinline ,(sprintf "  x%d = *p++;" (last lx)))))
	     (setq r (nconc1 
		      r 
		      `(cinline 
			,(concat 
			  (sprintf "*r++ += w%d * x%d " (car lw) (car lx))
			  (apply concat 
				 (all ((w (cdr lw)) (x (cdr lx)))
				   (sprintf " + w%d * x%d " w x)))
			  ";"))))
	     (setq lx (cdr (nconc1 lx (car lx)))))
	   r)
	 `((cinline "} end: ;")))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(de idx-m1fastconvol-check (kernel0 input0 output0 size)
  ((-int-) kernel0 input0 output0 size)
  (when (<> size kernel0)
    (error "kernel has wrong size"))
  (when (<> (- input0 output0) (1- size))
    (error "inconsistant sizes for convolution")) ())

(de error-notcontiguous ()
    (error "input or output not contiguous") ())

(dmd midx-m1fastconvol-check (input kernel output size)
  `(idx-m1fastconvol-check 
    (idx-dim ,kernel 0) (idx-dim ,input 0) (idx-dim ,output 0) ,size))

#? (idx-f1fastconvol2acc <input> <kernel> <output>) 
;; super-fast 1D convolution with size 2 kernel
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing floats.
(de idx-f1fastconvol2acc (input kernel output) 
    ((-idx1- (-float-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 2)
    (midx-m1fastconvolacc input kernel output 2 "float") ())

#? (idx-d1fastconvol2acc <input> <kernel> <output>) 
;; super-fast 1D convolution with size 2 kernel
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing doubles.
(de idx-d1fastconvol2acc (input kernel output) 
    ((-idx1- (-double-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 2)
    (midx-m1fastconvolacc input kernel output 2 "double") ())

#? (idx-u1fastconvol2acc <input> <kernel> <output>) 
;; super-fast 1D convolution with size 2 kernel
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing ubytes.
(de idx-u1fastconvol2acc (input kernel output) 
    ((-idx1- (-ubyte-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 2)
    (midx-m1fastconvolacc input kernel output 2 "char") ())


#? (idx-f1fastconvol3acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 3 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing floats.
(de idx-f1fastconvol3acc (input kernel output) 
    ((-idx1- (-float-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 3)
    (midx-m1fastconvolacc input kernel output 3 "float") ())

#? (idx-d1fastconvol3acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 3 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing doubles.
(de idx-d1fastconvol3acc (input kernel output) 
    ((-idx1- (-double-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 3)
    (midx-m1fastconvolacc input kernel output 3 "double") ())

#? (idx-u1fastconvol3acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 3 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing ubytes.
(de idx-u1fastconvol3acc (input kernel output) 
    ((-idx1- (-ubyte-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 3)
    (midx-m1fastconvolacc input kernel output 3 "char") ())



#? (idx-f1fastconvol4acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 4 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing floats.
(de idx-f1fastconvol4acc (input kernel output) 
    ((-idx1- (-float-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 4)
    (midx-m1fastconvolacc input kernel output 4 "float") ())

#? (idx-d1fastconvol4acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 4 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing doubles.
(de idx-d1fastconvol4acc (input kernel output) 
    ((-idx1- (-double-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 4)
    (midx-m1fastconvolacc input kernel output 4 "double") ())

#? (idx-u1fastconvol4acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 4 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing ubytes.
(de idx-u1fastconvol4acc (input kernel output) 
    ((-idx1- (-ubyte-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 4)
    (midx-m1fastconvolacc input kernel output 4 "char") ())


#? (idx-f1fastconvol5acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 5 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing floats.
(de idx-f1fastconvol5acc (input kernel output) 
    ((-idx1- (-float-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 5)
    (midx-m1fastconvolacc input kernel output 5 "float") ())

#? (idx-d1fastconvol5acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 5 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing doubles.
(de idx-d1fastconvol5acc (input kernel output) 
    ((-idx1- (-double-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 5)
    (midx-m1fastconvolacc input kernel output 5 "double") ())

#? (idx-u1fastconvol5acc <input> <kernel> <output>)
;; super-fast 1D convolution with size 5 kernel.
;; result is ACCUMULATED in the output
;; (with unrolled loops and optimal register assignment)
;; <input> and <output> must be contiguous idx1, <kernel> must be an idx1
;; all containing ubytes.
(de idx-u1fastconvol5acc (input kernel output) 
    ((-idx1- (-ubyte-)) input kernel output)
    (when (or (not (contiguousp input)) (not (contiguousp output)))
      (error-notcontiguous))
    (midx-m1fastconvol-check input kernel output 5)
    (midx-m1fastconvolacc input kernel output 5 "char") ())


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; THESE ARE THE "SLOW" CONVOLUTIONS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
#? *** Generic Convolutions
;; Regular 1D and 2D convolution and backconvolutions
;; for idx of doubles, floats, and ubytes.
;; Each function take 3 arguments: input, kernel, output.
;; The size of the output plus the size of the kernel
;; minus 1 must equal to the size of the input in
;; all dimensions (to prevent annoying boundary conditions).
;; The backconvolution can be seen the transposed operator
;; of the convolution.

#? (midx-m1convol <input> <kernel> <output>)
;; macro for 1D convolution. all arguments are idx1 (of any numerical type).
(dmd midx-m1convol (in kernel out)
  `(let ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1)))
     (idx-m2dotm1 uin ,kernel ,out)))

#? (idx-f1convol <input> <kernel> <output>)
;; 1D convolution on idx1 of floats
(de idx-f1convol (in kernel out) 
    ((-idx1- (-float-)) in kernel out)
    (midx-m1convol in kernel out))

#? (idx-m1convol <input> <kernel> <output>)
;; 1D convolution on idx1 of floats
;; (alias of idx-f1convol)
(setq idx-m1convol idx-f1convol)

#? (idx-d1convol <input> <kernel> <output>)
;; 1D convolution on idx1 of double
(de idx-d1convol (in kernel out) 
    ((-idx1- (-double-)) in kernel out)
    (midx-m1convol in kernel out))

#? (idx-u1convol <input> <kernel> <output>)
;; 1D convolution on idx1 of ubytes
(de idx-u1convol (in kernel out) 
    ((-idx1- (-ubyte-)) in kernel out)
    (midx-m1convol in kernel out))



#? (midx-m1backconvol <input> <kernel> <output>)
;; macro for 1D back convolution 
(dmd midx-m1backconvol (out kernel in)
  `(progn 
     (idx-m1clear ,in)
     (let ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1)))
       (idx-m1extm1acc ,out ,kernel uin))))

#? (idx-f1backconvol <input> <kernel> <output>)
;; 1D back convolution on floats
(de idx-f1backconvol (out kernel in) 
  ((-idx1- (-float-)) in kernel out)
  (midx-m1backconvol out kernel in))

#? (idx-m1backconvol <input> <kernel> <output>)
;; 1D back convolution on floats
;; (alias of idx-f1backconvol.
(setq idx-m1backconvol idx-f1backconvol)

#? (idx-d1backconvol <input> <kernel> <output>)
;; 1D back convolution on doubles
(de idx-d1backconvol (out kernel in) 
  ((-idx1- (-double-)) in kernel out)
  (midx-m1backconvol out kernel in))

#? (idx-u1backconvol <input> <kernel> <output>)
;; 1D back convolution on ubytes
(de idx-u1backconvol (out kernel in) 
  ((-idx1- (-ubyte-)) in kernel out)
  (midx-m1backconvol out kernel in))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
#? (midx-m2convol <input> <kernel> <output>)
;; macro for 2D convolution. all arguments are idx2.
(dmd midx-m2convol (in kernel out)
  `(let* ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1))
	  (uuin (unfold uin 1 (idx-dim ,kernel 1) 1)))
     (idx-m4dotm2 uuin ,kernel ,out)))

#? (idx-f2convol <input> <kernel> <output>)
;; 2D convolution. all arguments are idx2 of floats
(de idx-f2convol (in kernel out) 
  ((-idx2- (-float-)) in kernel out)
  (midx-m2convol in kernel out))

#? (idx-m2convol <input> <kernel> <output>)
;; 2D convolution. all arguments are idx2 of floats.
;; (alias for idx-f2convol)
(setq idx-m2convol idx-f2convol)

#? (idx-d2convol <input> <kernel> <output>)
;; 2D convolution. all arguments are idx2 of doubles.
(de idx-d2convol (in kernel out) 
  ((-idx2- (-double-)) in kernel out)
  (midx-m2convol in kernel out))

#? (idx-u2convol <input> <kernel> <output>)
;; 2D convolution. all arguments are idx2 of ubytes
(de idx-u2convol (in kernel out) 
  ((-idx2- (-ubyte-)) in kernel out)
  (midx-m2convol in kernel out))


#? (midx-m2convolacc <input> <kernel> <output>)
;; macro for 2D convolution with accumulation . all arguments are idx2.
(dmd midx-m2convolacc (in kernel out)
  `(let* ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1))
	  (uuin (unfold uin 1 (idx-dim ,kernel 1) 1)))
     (idx-m4dotm2acc uuin ,kernel ,out)))

#? (idx-f2convolacc <input> <kernel> <output>)
;; 2D convolution with accumulation . all arguments are idx2 of floats
(de idx-f2convolacc (in kernel out) 
      ((-idx2- (-float-)) in kernel out)
      (midx-m2convolacc in kernel out))

#? (idx-m2convolacc <input> <kernel> <output>)
;; 2D convolution with accumulation . all arguments are idx2.
(setq idx-m2convolacc idx-f2convolacc)

#? (idx-d2convolacc <input> <kernel> <output>)
;; 2D convolution with accumulation . all arguments are idx2 of double.
(de idx-d2convolacc (in kernel out) 
      ((-idx2- (-double-)) in kernel out)
      (midx-m2convolacc in kernel out))

#? (idx-u2convolacc <input> <kernel> <output>)
;; 2D convolution with accumulation . all arguments are idx2 of ubyte.
(de idx-u2convolacc (in kernel out) 
      ((-idx2- (-ubyte-)) in kernel out)
      (midx-m2convolacc in kernel out))


#? (midx-m2backconvol <input> <kernel> <output>)
;; macro for 2D back-convolution. all arguments are idx2.
(dmd midx-m2backconvol (out kernel in)
  `(progn 
     (idx-m2clear ,in)
     (let* ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1))
	    (uuin (unfold uin 1 (idx-dim ,kernel 1) 1)))
       (idx-m2extm2acc ,out ,kernel uuin))))

#? (idx-f2backconvol <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-f2backconvol (out kernel in) 
  ((-idx2- (-float-)) in kernel out)
  (midx-m2backconvol out kernel in))

#? (idx-m2backconvol <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats
;; (alias for idx-f2baconvol)
(setq idx-m2backconvol idx-f2backconvol)

#? (idx-d2backconvol <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-d2backconvol (out kernel in) 
  ((-idx2- (-double-)) in kernel out)
  (midx-m2backconvol out kernel in))

#? (idx-u2backconvol <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-u2backconvol (out kernel in) 
  ((-idx2- (-ubyte-)) in kernel out)
  (midx-m2backconvol out kernel in))


#? (midx-m2backconvolacc <input> <kernel> <output>)
;; macro for 2D back-convolution. all arguments are idx2.
(dmd midx-m2backconvolacc (out kernel in)
  `(progn 
     (let* ((uin (unfold ,in 0 (idx-dim ,kernel 0) 1))
	    (uuin (unfold uin 1 (idx-dim ,kernel 1) 1)))
       (idx-m2extm2acc ,out ,kernel uuin))))

#? (idx-f2backconvolacc <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-f2backconvolacc (out kernel in) 
  ((-idx2- (-float-)) in kernel out)
  (midx-m2backconvolacc out kernel in))

#? (idx-m2backconvolacc <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(setq idx-m2backconvolacc idx-f2backconvolacc)

#? (idx-d2backconvolacc <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-d2backconvolacc (out kernel in) 
  ((-idx2- (-double-)) in kernel out)
  (midx-m2backconvolacc out kernel in))

#? (idx-u2backconvolacc <input> <kernel> <output>)
;; 2D back-convolution. all arguments are idx2 of floats.
(de idx-u2backconvolacc (out kernel in) 
  ((-idx2- (-ubyte-)) in kernel out)
  (midx-m2backconvolacc out kernel in))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? *** 1D and 2D Subsampling and Oversampling
;; functions that subsample 1D and 2D idx
;; where each output sample is the average of <n>
;; adjacent input samples (with no overlap).
;; The oversmapling functions simply replicate
;; the samples <n> times.

#? (midx-m1subsample <input> <n> <output>)
;; macro for subsampling and summing of size <n>
(dmd midx-m1subsample (in n out)
  `(let ((uin (unfold ,in 0 ,n ,n)))
     (idx-m1clear ,out)
     (idx-eloop ((z uin))
       (idx-m1add z ,out ,out))))

#? (idx-f1subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-f1subsample (in n out) 
  ((-int-) n) 
  ((-idx1- (-float-)) in out)
  (midx-m1subsample in n out))

#? (idx-m1subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(setq idx-m1subsample idx-f1subsample)

#? (idx-d1subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-d1subsample (in n out) 
  ((-int-) n) 
  ((-idx1- (-double-)) in out)
  (midx-m1subsample in n out))

#? (idx-u1subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-u1subsample (in n out) 
  ((-int-) n) 
  ((-idx1- (-ubyte-)) in out)
  (midx-m1subsample in n out))



#? (midx-m1oversample <input> <n> <output>)
;; macro for oversampling of size <n>
(dmd midx-m1oversample (out n in)
  `(let ((uin (unfold ,in 0 ,n ,n)))
     (idx-eloop ((z uin))
       (idx-copy ,out z)) ()))

#? (idx-f1oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-f1oversample (out n in) 
  ((-int-) n) 
  ((-idx1- (-float-)) in out)
  (midx-m1oversample out n in))

#? (idx-m1oversample <input> <n> <output>)
;; oversampling of size <n>
(setq idx-m1oversample idx-f1oversample) 

#? (idx-d1oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-d1oversample (out n in) 
  ((-int-) n) 
  ((-idx1- (-double-)) in out)
  (midx-m1oversample out n in))

#? (idx-u1oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-u1oversample (out n in) 
  ((-int-) n) 
  ((-idx1- (-ubyte-)) in out)
  (midx-m1oversample out n in))


#? (midx-m2subsample <input> <n> <output>)
;; macro for subsampling and summing of size <n>
(dmd midx-m2subsample (in nlin ncol out)
  `(let* ((uin (unfold ,in 0 ,nlin ,nlin))
	  (uuin (unfold uin 1 ,ncol ,ncol)))
     (idx-clear ,out)
     (idx-eloop ((z1 uuin))
       (idx-eloop ((z2 z1))
	 (idx-add z2 ,out ,out))) ()))

#? (idx-f2subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-f2subsample (in nlin ncol out) 
  ((-int-) nlin ncol)
  ((-idx2- (-float-)) in out)
  (midx-m2subsample in nlin ncol out))

#? (idx-m2subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(setq idx-m2subsample idx-f2subsample)

#? (idx-d2subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-d2subsample (in nlin ncol out) 
  ((-int-) nlin ncol)
  ((-idx2- (-double-)) in out)
  (midx-m2subsample in nlin ncol out))

#? (idx-u2subsample <input> <n> <output>)
;; subsampling and summing of size <n>
(de idx-u2subsample (in nlin ncol out) 
  ((-int-) nlin ncol)
  ((-idx2- (-ubyte-)) in out)
  (midx-m2subsample in nlin ncol out))



#? (midx-m2oversample <input> <n> <output>)
;; macro for oversampling of size <n>
(dmd midx-m2oversample (out nlin ncol in)
  `(let* ((uin (unfold ,in 0 ,nlin ,nlin))
	  (uuin (unfold uin 1 ,ncol ,ncol)))
     (idx-eloop ((z1 uuin))
       (idx-eloop ((z2 z1))
	 (idx-copy ,out z2))) ()))

#? (idx-f2oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-f2oversample (out nlin ncol in) 
  ((-int-) nlin ncol)
  ((-idx2- (-float-)) in out)
  (midx-m2oversample out nlin ncol in))

#? (idx-m2oversample <input> <n> <output>)
;; oversampling of size <n>
(setq idx-m2oversample idx-f2oversample)

#? (idx-d2oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-d2oversample (out nlin ncol in) 
  ((-int-) nlin ncol)
  ((-idx2- (-double-)) in out)
  (midx-m2oversample out nlin ncol in))

#? (idx-u2oversample <input> <n> <output>)
;; oversampling of size <n>
(de idx-u2oversample (out nlin ncol in) 
  ((-int-) nlin ncol)
  ((-idx2- (-ubyte-)) in out)
  (midx-m2oversample out nlin ncol in))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(dhc-make ()
	  error-notcontiguous
	  idx-m1fastconvol-check 
	  idx-f1fastconvol2acc 
	  idx-d1fastconvol2acc 
	  idx-u1fastconvol2acc 
	  idx-f1fastconvol3acc 
	  idx-d1fastconvol3acc 
	  idx-u1fastconvol3acc 
	  idx-f1fastconvol4acc 
	  idx-d1fastconvol4acc 
	  idx-u1fastconvol4acc 
	  idx-f1fastconvol5acc 
	  idx-d1fastconvol5acc 
	  idx-u1fastconvol5acc 
	  idx-f1convol 
	  idx-d1convol 
	  idx-u1convol 
	  idx-f1backconvol 
	  idx-d1backconvol 
	  idx-u1backconvol 
	  idx-f2convol 
	  idx-d2convol 
	  idx-u2convol 
	  idx-f2convolacc 
	  idx-d2convolacc 
	  idx-u2convolacc 
	  idx-f2backconvol 
	  idx-d2backconvol 
	  idx-u2backconvol 
	  idx-f2backconvolacc 
	  idx-d2backconvolacc 
	  idx-u2backconvolacc 
	  idx-f1subsample 
	  idx-d1subsample 
	  idx-u1subsample 
	  idx-f1oversample 
	  idx-d1oversample 
	  idx-u1oversample 
	  idx-f2subsample 
	  idx-d2subsample 
	  idx-u2subsample 
	  idx-f2oversample 
	  idx-d2oversample 
	  idx-u2oversample 
	  )