idx.hpp

Gaussian Mixture Algorithm

/***************************************************************************
 *   Copyright (C) 2008 by Yann LeCun   *
 *   [e-mail]   *
 *   All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Redistribution under a license not approved by the Open Source
 *       Initiative (http://www.opensource.org) must display the
 *       following acknowledgement in all advertising material:
 *        This product includes software developed at the Courant
 *        Institute of Mathematical Sciences (http://cims.nyu.edu).
 *     * The names of the authors may not be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL ThE AUTHORS BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ***************************************************************************/

#ifndef Idx_HPP
#define Idx_HPP

#include <sstream>
#include <vector>

namespace ebl {

////////////////////////////////////////////////////////////////

//! Idx elements and dimensions error checking macros

//! Calls ylerror if src0 and src1 have different number of elements.
#define idx_checknelems2_all(src0, src1) \
  if ((src0).nelements() != (src1).nelements()) ylerror("incompatible Idx sizes\n");

//! Calls ylerror if src0 and src1 and src2 have different number of elements.
#define idx_checknelems3_all(src0, src1, src2) \
  if (((src0).nelements() != (src1).nelements()) || \
  		((src0).nelements() != (src2).nelements())) \
  		ylerror("incompatible Idx sizes\n");

//! Calls ylerror if src0 and o0 do not match.
#define idx_checkorder1(src0, o0) \
  if ((src0).order() != o0) \
    ylerror("Idx has wrong order");

//! Calls ylerror if src0,src1 and o0,o1 do not match.
#define idx_checkorder2(src0, o0, src1, o1) \
  if (((src0).order() != o0) || ((src1).order() != o1)) \
    ylerror("Idx have incompatible orders");

//! Calls ylerror if src0,src1,src2 and o0,o1,o2 do not match.
#define idx_checkorder3(src0, o0, src1, o1, src2, o2) \
  if (((src0).order() != o0) || ((src1).order() != o1) || ((src2).order() != o2)) \
    ylerror("Idx have incompatible orders");

//! Calls ylerror if src0.order(), src1.order() and src2.order() differ
#define idx_checkorder3_all(src0, src1, src2) \
  if (((src0).order() != (src1).order()) || ((src0).order() != (src2).order())) \
    ylerror("Idx have incompatible orders");

//! Calls ylerror if src0.dim(0) and src1.dim(d1) don't match e0,e1
#define idx_checkdim2(src0, d0, e0, src1, d1, e1) \
	if (((src0).dim(d0) != e0) || ((src1).dim(d1) != e1))\
    ylerror("Idx have incompatible dimensions");

//! Calls ylerror if src0.dim(d) and src1.dim(d) don't match
#define idx_checkdim2_all(src0, src1, d) \
	if ((src0).dim(d) != (src1).dim(d)) \
    ylerror("Idx have incompatible dimensions");

//! Calls ylerror if src0.dim(d) and src1.dim(d) and src2.dim(d) don't match
#define idx_checkdim3_all(src0, src1, src2, d) \
	if (((src0).dim(d) != (src1).dim(d)) || \
			((src0).dim(d) != (src2).dim(d))) \
    ylerror("Idx have incompatible dimensions");

//! Calls ylerror if src0.dim(d) and src1.dim(d) and src2.dim(d)
//! and src3.dim(d) don't match
#define idx_checkdim4_all(src0, src1, src2, src3, d) \
	if (((src0).dim(d) != (src1).dim(d)) || \
			((src0).dim(d) != (src2).dim(d)) || \
			((src0).dim(d) != (src3).dim(d))) \
    ylerror("Idx have incompatible dimensions");

////////////////////////////////////////////////////////////////

// TODO: these macros are untested (YLC)

//! cidxN_bloopX: macros to loop simultaneously over elements 
//! of X Idx'es of order at least N. Can be used as follows:
//! { double *z0, *z1; 
//!   intg i;
//!   cidx1_bloop2(i, z0, myidx0, z1, myidx1) { *z0 = *z1 + 4; }
//! }
//! { float *z0;
//!   intg i,j;
//!   cidx2_bloop1(i, j, z0, myidx0) { *z0 *= 2; }
//! }
//! Variable i is a loop index, myidx0 is an Idx of any type whose
//! order must be at least 1, and and z0 is a pointer to the numerical 
//! type of myidx0.
//! It is best to enclose each cidx1_bloopX in its own brace scope
//! because the macro creates temprary variables with a fixed name.
//! These macros should be used over the idx_aloopX macros whenever
//! possible, because they are considerably more efficient.
//! unlike the aloop macros, these macros manipulate regular pointers 
//! with simple incrementation, as opposed to iterators with complicated 
//! logic.
#define cidx1_bloop1(i,p0,src0)		\
  if ((src0).order() < 1) ylerror("Idx has wrong order"); \
  intg _n0 = (src0).dim(0), _m0 = (src0).mod(0);		\
  for (i=0, p0=(src0).idx_ptr(); i<_n0; i++, p0+=_m0) 

#define cidx1_bloop2(i,p0,src0,p1,src1)		 \
  if (((src0).order() < 1)||((src1).order() < 1))  ylerror("Idx has wrong order"); \
  intg _n0 = (src0).dim(0), _m0 = (src0).mod(0); _m1 = (src1).mod(0);	 \
  idx_checkdim2_all(src0,src1,0)					\
  for (i=0, p0=(src0).idx_ptr(), p1=(src1).idx_ptr(); i<_n0; i++, p0+=_m0, p1+=_m1) 

#define cidx1_bloop3(i,p0,src0,p1,src1,p2,src2)		 \
  intg _n0 = (src0).dim(0), _m0 = (src0).mod(0), _m1 = (src1).mod(0), _m2 = (src2).mod(0); \
  idx_checkdim3_all(src0,src1,src2,0)					\
  for (i=0, p0=(src0).idx_ptr(), p1=(src1).idx_ptr(), p2=(src2).idx_ptr(); i<_n0; i++, p0+=_m0, p1+=_m1, p2+=_m2) 

#define cidx1_bloop4(i,p0,src0,p1,src1,p2,src2,p3,src3)	 \
  intg _n0 = (src0).dim(0), _m0 = (src0).mod(0), _m1 = (src1).mod(0), _m2 = (src2).mod(0), _m3 = (src3).mod(0); \
    idx_checkdim4_all(src0,src1,src2,src3,0)				\
    for (i=0, p0=(src0).idx_ptr(), p1=(src1).idx_ptr(), p2=(src2).idx_ptr(), p3=(src3).idx_ptr(); i<_n0; i++, p0+=_m0, p1+=_m1, p2+=_m2, p3+=_m3) 


#define cidx2_bloop1(i,j,p0,src0)			  \
  if ((src0).order() < 2) ylerror("Idx has wrong order"); \
  intg _n00 = (src0).dim(0), _m00 = (src0).mod(0);		\
  intg _n01 = (src0).dim(1), _m01 = (src0).mod(1);		\
  for (i=0, p0=(src0).idx_ptr(); i<_n00; i++, p0+=_m00-_n01*_m01) \
  for (j=0; i<_n01; j++, p0+=_m01) 

#define cidx2_bloop2(i,j,p0,src0,p1,src1)		  \
  if ((src0).order() < 2) ylerror("Idx has wrong order"); \
  intg _n00 = (src0).dim(0), _m00 = (src0).mod(0);		\
  intg _n01 = (src0).dim(1), _m01 = (src0).mod(1);		\
  intg _n10 = (src1).dim(0), _m10 = (src1).mod(0);		\
  intg _n11 = (src1).dim(1), _m11 = (src1).mod(1);		\
  idx_checkdim2_all(src0,src1,0)					\
  idx_checkdim2_all(src0,src1,1)					\
  for (i=0, p0=(src0).idx_ptr(), p1=(src1).idx_ptr(); i<_n00; i++, p0+=_m00-_n01*_m01, p1+=_m10-_n11*_m11) \
    for (j=0; i<_n01; j++, p0+=_m01, p1+=_m11) 

#define cidx2_bloop3(i,j,p0,src0,p1,src1,p2,src2)		  \
  if ((src0).order() < 2) ylerror("Idx has wrong order"); \
  intg _n00 = (src0).dim(0), _m00 = (src0).mod(0);		\
  intg _n01 = (src0).dim(1), _m01 = (src0).mod(1);		\
  intg _n10 = (src1).dim(0), _m10 = (src1).mod(0);		\
  intg _n11 = (src1).dim(1), _m11 = (src1).mod(1);		\
  intg _n20 = (src2).dim(0), _m20 = (src2).mod(0);		\
  intg _n21 = (src2).dim(1), _m21 = (src2).mod(1);		\
  idx_checkdim3_all(src0,src1,src2,0)					\
  idx_checkdim3_all(src0,src1,src2,1)					\
  for (i=0, p0=(src0).idx_ptr(), p1=(src1).idx_ptr(), p2=(src2).idx_ptr(); i<_n00; i++, p0+=_m00-_n01*_m01, p1+=_m10-_n11*_m11, p2+=_m20-_n21*_m21) \
    for (j=0; i<_n01; j++, p0+=_m01, p1+=_m11, p2+=_m21) 


////////////////////////////////////////////////////////////////

//! call these macros like this:
//! { idx_bloop1(la, a) { r += la.get(); } }
//! exmaple: matrix-vector product  a x b -> c
//! { idx_bloop3(la, a, lb, b, lc, c) { dot_product(la,lb,lc); } }
//! It's advisable to encase every bloop in its own scope to
//! prevent name clashes if the same loop variable is used
//! multiple times in a scope.

// bloop/eloop macros

// Okay, either C++ really suxx0rz or I'm a really st00pid Lisp-head.
// Why can't we define syntax-changing macros like in Lisp? Basically,
// I can't encapsulate the allocation
// of temporary variable for a loop inside the macro, unless
// I define a "begin" macro and an "end" macro. I would like
// to define idx_bloop so I can do:
// idx_bloop2(lm, m, lv, v) { idx_dot(lm,lm,v); }
// but I can't do that because I have to allocate lm and lv
// inside the macro, hence I need to know the type.
// Now the call would be:
// idx_bloop2(lm, m, double, lv, v, double) { idx_dot(lm,lm,v); }
// But that still doesn't quite work because if I declare lm and lv
// then I can't reuse the same symbols for another loop in the same
// scope. The only way out is to force the user to encase every
// bloop call inside braces, or to not reuse the same synbol twice
// for a looping Idx. I thought about generating a mangled name
// but couldn't find a way to make it useful.
// If a macro could define its own scope that would be great.


#if USING_STL_ITERS == 1

#define idx_bloop1(dst0,src0,type0)	 \
  for ( DimIter<type0> dst0(src0,0); dst0.notdone(); ++dst0)

#define idx_bloop2(dst0,src0,type0,dst1,src1,type1)				\
  idx_checkdim2_all(src0, src1, 0); \
  DimIter<type0> dst0(src0,0); \
  DimIter<type1> dst1(src1,0); \
  for ( ; dst0.notdone(); ++dst0, ++dst1)

#define idx_bloop3(dst0,src0,type0,dst1,src1,type1,dst2,src2,type2)	\
  idx_checkdim3_all(src0, src1, src2, 0); \
  DimIter<type0> dst0(src0,0);					\
  DimIter<type1> dst1(src1,0);					\
  DimIter<type2> dst2(src2,0);					\
  for ( ; dst0.notdone(); ++dst0, ++dst1, ++dst2)

#define idx_bloop4(dst0,src0,type0,dst1,src1,type1,dst2,src2,type2,dst3,src3,type3) \
  idx_checkdim4_all(src0, src1, src2, src3, 0); \
  DimIter<type0> dst0(src0,0);					\
  DimIter<type1> dst1(src1,0);					\
  DimIter<type2> dst2(src2,0);					\
  DimIter<type3> dst3(src3,0);					\
  for ( ; dst0.notdone(); ++dst0, ++dst1, ++dst2, ++dst3)

// eloop macros

#define idx_eloop1(dst0,src0,type0) \
  DimIter<type0> dst0(src0,src0.order()-1);	\
  for ( ; dst0.notdone(); ++dst0)

#define idx_eloop2(dst0,src0,type0,dst1,src1,type1)				\
  if ( (src0).dim((src0).order()) != (src1).dim((src1).order()) ) ylerror("incompatible Idxs for bloop\n"); \
  DimIter<type0> dst0(src0,(src0).order()-1); \
  DimIter<type1> dst1(src1,(src1).order()-1); \
  for ( ; dst0.notdone(); ++dst0, ++dst1)


////////////////////////////////////////////////////////////////
// aloop macros: loop over all elements of an Idx
// These macros are somewhat inefficient and should be used as little
// as possible, or whenever simplicity is preferable to speed.

// Loops over all elements of an idx. This takes a pointer to
// the data type of idx elements, and a blank IdxIter object:
// idx_aloop1(data_pointer,idxiter,&idx) { do_stuff(data_pointer); }
// Example of use: add 1 to all element of m:
//  Idx<double> m(3,4);
//  ScalarIter<double> p;
//  idx_aloop1(p,&m) { *p += 1; }
#define idx_aloop1_on(itr0,src0)					\
  for ( ; itr0.notdone(); ++itr0)


// this loops simultaneously over all elements of 2 Idxs.
// The two Idxs can have different structures as long as they have
// the same total number of elements.
#define idx_aloop2_on(itr0,src0,itr1,src1)			\
	idx_checknelems2_all(src0, src1); \
  for ( ; itr0.notdone(); ++itr0, ++itr1)


#define idx_aloop3_on(itr0,src0,itr1,src1,itr2,src2)	\
	idx_checknelems3_all(src0, src1, src2); \
  for (; itr0.notdone(); ++itr0, ++itr1, ++itr2)


// high level aloop macros.
// These should be enclosed in braces, to avoid name clashes
#define idx_aloop1(itr0,src0,type0) \
  ScalarIter<type0> itr0(src0); \
  idx_aloop1_on(itr0,src0)


#define idx_aloop2(itr0,src0,type0,itr1,src1,type1)		\
  ScalarIter<type0> itr0(src0); \
  ScalarIter<type1> itr1(src1); \
  idx_checknelems2_all(src0, src1); \
  for (; itr0.notdone(); ++itr0, ++itr1)

#define idx_aloop3(itr0,src0,type0,itr1,src1,type1,itr2,src2,type2) \
  ScalarIter<type0> itr0(src0); \
  ScalarIter<type1> itr1(src1); \
  ScalarIter<type2> itr2(src2); \
  idx_checknelems3_all(src0, src1, src2); \
  for (; itr0.notdone(); ++itr0, ++itr1, ++itr2)

#else

#define idx_bloop1(dst0,src0,type0)	 \
  IdxLooper<type0> dst0(src0,0); \
  for ( ; dst0.notdone(); dst0.next())

#define idx_bloop2(dst0,src0,type0,dst1,src1,type1)				\
  idx_checkdim2_all(src0, src1, 0); \
  IdxLooper<type0> dst0(src0,0); \
  IdxLooper<type1> dst1(src1,0); \
  for ( ; dst0.notdone(); dst0.next(), dst1.next())

#define idx_bloop3(dst0,src0,type0,dst1,src1,type1,dst2,src2,type2)	\
  idx_checkdim3_all(src0, src1, src2, 0); \
  IdxLooper<type0> dst0(src0,0);					\
  IdxLooper<type1> dst1(src1,0);					\
  IdxLooper<type2> dst2(src2,0);					\
  for ( ; dst0.notdone(); dst0.next(), dst1.next(), dst2.next())