idx.hpp

Gaussian Mixture Algorithm

//template <typename T>
//Idx<T>::Idx( const Idx<T>& other )
//	:storage(other.storage),
//	 spec(other.spec)
// {
//	storage->lock();
// }

template <class T> Idx<T>::Idx(Srg<T> *srg, IdxSpec &s) {
  spec = s;
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(Srg<T> *srg, intg o) : spec(o) {
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx() : spec(0) {
  storage = new Srg<T>();
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(Srg<T> *srg, intg o, intg size0) : spec(o,size0) {
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(intg size0) : spec(0,size0) {
  storage = new Srg<T>();
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(Srg<T> *srg, intg o, intg size0, intg size1) : spec(o,size0,size1) {
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(intg size0, intg size1) : spec(0,size0,size1) {
  storage = new Srg<T>();
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(Srg<T> *srg, intg o, intg size0, intg size1, intg size2) : spec(o,size0,size1,size2) {
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(intg size0, intg size1, intg size2) : spec(0,size0,size1,size2) {
  storage = new Srg<T>();
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(Srg<T> *srg, intg o, intg s0, intg s1, intg s2, intg s3, intg s4, intg s5, intg s6, intg s7) : spec(o,s0,s1,s2,s3,s4,s5,s6,s7) {
  storage = srg;
  growstorage();
  storage->lock();
}

template <class T> Idx<T>::Idx(intg s0, intg s1, intg s2, intg s3, intg s4, intg s5, intg s6, intg s7) : spec(0,s0,s1,s2,s3,s4,s5,s6,s7) {
  storage = new Srg<T>();
  growstorage();
  storage->lock();
}


template <class T> intg Idx<T>::setoffset(intg o) {
  if (o<0) { ylerror("Idx::setoffset: offset must be positive"); }
  if (o > spec.offset) {
    spec.setoffset(o);
    growstorage();
    return o;
  } else {
    spec.setoffset(o);
    return o;
  }
}

template <class T> void Idx<T>::resize(intg s0, intg s1, intg s2, intg s3, intg s4, intg s5, intg s6, intg s7) {
  spec.resize(s0,s1,s2,s3,s4,s5,s6,s7);
  growstorage();
}

template <class T> void Idx<T>::resize_chunk(intg s_chunk, intg s0, intg s1, intg s2, intg s3, intg s4, intg s5, intg s6, intg s7) {
  spec.resize(s0,s1,s2,s3,s4,s5,s6,s7);
  growstorage_chunk(s_chunk);
}

//template<typename T, typename SizeIter>
//void Idx<T>::resize( SizeIter& sizesBegin, SizeIter& sizesEnd ){
//	const int ndims = std::distance( sizesBegin, sizesEnd );
//	if ( ndims > MAXDIMS ){
//		std::ostringstream oss;
//		oss<<"Number of dimensions ("<<ndims<<") exceeds MAX_DIMS ("<<MAXDIMS<<")";
//		ylerror(oss.str());
//	}
//
//	std::vector<T> sizes(ndims+1);
//	std::copy(sizesBegin, sizesEnd, sizes.begin());
//	sizes.back() = -1;
//	spec.resize( sizesBegin, sizesEnd );
//	growstorage();
//}


template <class T> Idx<T> Idx<T>::select(int d, intg i) {
  Idx<T> r(storage,spec.getoffset());
  spec.select_into(&r.spec, d, i);
  return r;
}

template <class T> Idx<T> Idx<T>::narrow(int d, intg s, intg o) {
  Idx<T> r(storage,spec.getoffset());
  spec.narrow_into(&r.spec, d, s, o);
  return r;
}

template <class T> Idx<T> Idx<T>::transpose(int d1, int d2) {
  Idx<T> r(storage,spec.getoffset());
  spec.transpose_into(&r.spec, d1, d2);
  return r;
}

template <class T> Idx<T> Idx<T>::transpose(int *p) {
  Idx<T> r(storage,spec.getoffset());
  spec.transpose_into(&r.spec, p);
  return r;
}

template <class T> Idx<T> Idx<T>::unfold(int d, intg k, intg s) {
  Idx<T> r(storage,spec.getoffset());
  spec.unfold_into(&r.spec, d, k, s);
  return r;
}


////////////////////////////////////////////////////////////////
// access methods

// get element of Idx1
template <class T> T *Idx<T>::ptr(intg i0) {
  if (spec.ndim != 1) ylerror("not an Idx1");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  return storage->data + spec.offset + i0*spec.mod[0];
}

// get element of Idx2
template <class T> T *Idx<T>::ptr(intg i0, intg i1) {
  if (spec.ndim != 2) ylerror("not an Idx2");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  return storage->data + spec.offset + i0*spec.mod[0] + i1*spec.mod[1];
}

// get element of Idx3
template <class T> T *Idx<T>::ptr(intg i0, intg i1, intg i2) {
  if (spec.ndim != 3) ylerror("not an Idx3");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  if ((i2 < 0) || (i2 >= spec.dim[2])) ylerror("index 2 out of bound");
  return storage->data + spec.offset + i0*spec.mod[0] + i1*spec.mod[1] + i2*spec.mod[2];
}


// return a pointer to an element of an Idx
// generic function for order>3
template <class T> T *Idx<T>::ptr(intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7) {
  try {
    // check that we passed the right number of indices
    // and that they are all positive
    switch (spec.ndim) {
    case 8: if (i7 < 0) throw(-8);break;
    case 7: if ((i6 < 0) || (i7 != -1)) throw(-7);break;
    case 6: if ((i5 < 0) || (i6 != -1)) throw(-6);break;
    case 5: if ((i4 < 0) || (i5 != -1)) throw(-5);break;
    case 4: if ((i3<0)||(i2<0)||(i1<0)||(i0<0)||(i4 != -1)) throw(-4);break;
    default: throw(10);
    }
    // now compute offset, and check that all
    // indices are within bounds.
    intg k = 0;
    switch (spec.ndim) {
    case 8: k += spec.mod[7]*i7; if (i7 >= spec.dim[7])  throw(7);
    case 7: k += spec.mod[6]*i6; if (i6 >= spec.dim[6])  throw(6);
    case 6: k += spec.mod[5]*i5; if (i5 >= spec.dim[5])  throw(5);
    case 5: k += spec.mod[4]*i4; if (i4 >= spec.dim[4])  throw(4);
    case 4: k += spec.mod[3]*i3; if (i3 >= spec.dim[3])  throw(3);
    }
    k += spec.mod[2]*i2; if (i2 >= spec.dim[2])  throw(2);
    k += spec.mod[1]*i1; if (i1 >= spec.dim[1])  throw(1);
    k += spec.mod[0]*i0; if (i0 >= spec.dim[0])  throw(0);
    return storage->data + spec.offset + k;
  }
  catch(int k) {
    if (k==10) ylerror("Idx::get: number of indices and order are different");
    if (k < 0) {
      ylerror("Idx::get: wrong number of indices, or negative index");
    } else {
      ylerror("Idx::get: index out of bound");
    }
    return NULL;
  }
}

////////////////////////////////////////////////////////////////
// get

// get element of Idx0
template <class T> T Idx<T>::get() {
  if (spec.ndim != 0) ylerror("not an Idx0");
  return (storage->data)[spec.offset];
}

// get element of Idx1
template <class T> T Idx<T>::get(intg i0) {
  if (spec.ndim != 1) ylerror("not an Idx1");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0]];
}

// get element of Idx2
template <class T> T Idx<T>::get(intg i0, intg i1) {
  if (spec.ndim != 2) ylerror("not an Idx2");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0] + i1*spec.mod[1]];
}

// get element of Idx3
template <class T> T Idx<T>::get(intg i0, intg i1, intg i2) {
  if (spec.ndim != 3) ylerror("not an Idx3");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  if ((i2 < 0) || (i2 >= spec.dim[2])) ylerror("index 2 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0] + i1*spec.mod[1] + i2*spec.mod[2]];
}

// get element of an Idx of any order
template <class T> T Idx<T>::get(intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7) {
  return *ptr(i0,i1,i2,i3,i4,i5,i6,i7);
}

////////////////////////////////////////////////////////////////
// set

// set the element of Idx0
template <class T> T Idx<T>::set(T val) {
  if (spec.ndim != 0) ylerror("not an Idx0");
  return (storage->data)[spec.offset] = val;
}

// set the element of Idx1
template <class T> T Idx<T>::set(T val, intg i0) {
  if (spec.ndim != 1) ylerror("not an Idx1");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0]] = val;
}

// set the element of Idx2
template <class T> T Idx<T>::set(T val, intg i0, intg i1) {
  if (spec.ndim != 2) ylerror("not an Idx2");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0] + i1*spec.mod[1]] = val;
}

// set the element of Idx3
template <class T> T Idx<T>::set(T val, intg i0, intg i1, intg i2) {
  if (spec.ndim != 3) ylerror("not an Idx3");
  if ((i0 < 0) || (i0 >= spec.dim[0])) ylerror("index 0 out of bound");
  if ((i1 < 0) || (i1 >= spec.dim[1])) ylerror("index 1 out of bound");
  if ((i2 < 0) || (i2 >= spec.dim[2])) ylerror("index 2 out of bound");
  return (storage->data)[spec.offset + i0*spec.mod[0] + i1*spec.mod[1] + i2*spec.mod[2]] = val;
}


// set an element of an Idx of any order.
template <class T> T Idx<T>::set(T val, intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7) {
    return *ptr(i0,i1,i2,i3,i4,i5,i6,i7) = val;
}

template <class T> int Idx<T>::fdump(FILE *f) {
  if (spec.ndim == 0) {
    fprintf(f,"[@ %g]\n",this->get());
  } else if (spec.ndim == 1) {
    fprintf(f,"[");
    for (intg i=0; i<dim(0); i += mod(0)) {
      fprintf(f,"%g ",(storage->data)[spec.offset + i]);
    }
    fprintf(f,"]\n");
  } else {
    fprintf(f,"[");
    { idx_bloop1(p,*this,T) { p.fdump(f); } }
    fprintf(f,"]\n");
  }
  return 0;
}

////////////////////////////////////////////////////////////////
// STL-style iterator creators
template <typename T> typename Idx<T>::scalar_iterator Idx<T>::scalars_begin(){
	return scalar_iterator(*this);
}

template <typename T> typename Idx<T>::scalar_iterator Idx<T>::scalars_end(){
	return scalar_iterator(*this, false);
}

template <typename T> typename Idx<T>::reverse_scalar_iterator Idx<T>::scalars_rbegin(){
	return reverse_scalar_iterator(*this);
}

template <typename T> typename Idx<T>::reverse_scalar_iterator Idx<T>::scalars_rend(){
	return reverse_scalar_iterator(*this);
}

template <typename T>
typename Idx<T>::dimension_iterator
Idx<T>::dim_begin(int dd){
	return dimension_iterator(*this,dd);
}

template <typename T>
typename Idx<T>::dimension_iterator
Idx<T>::dim_end(int dd){
	return dimension_iterator(*this,dd,false);
}

template <typename T>
typename Idx<T>::reverse_dimension_iterator
Idx<T>::dim_rbegin(int dd){
	return reverse_dimension_iterator(*this,dd);
}

template <typename T>
typename Idx<T>::reverse_dimension_iterator
Idx<T>::dim_rend(int dd){
	return reverse_dimension_iterator(*this,dd,false);
}


#if USING_STL_ITERS == 0
////////////////////////////////////////////////////////////////
// an IdxLooper is an iterator for Idxs.
// It is actually a subclass of Idx.
// These are not C++ iterators in the classical sense.

template <class T> IdxLooper<T>::IdxLooper(Idx<T> &idx, int ld) : Idx<T>((dummyt*)0) {
  i = 0;
  dimd = idx.spec.dim[ld];
  modd = idx.spec.mod[ld];
  idx.spec.select_into(&(this->spec), ld, i);
  this->storage = idx.storage;
  this->storage->lock();
}

// like ++
// CAUTION: this doesn't do array bound checking
// because we coudn't use a for loop if it did.
template <class T> T *IdxLooper<T>::next() {
  i++;
  this->spec.offset += modd;
  return this->storage->data + this->spec.offset;
}

// return true when done.
template <class T> bool IdxLooper<T>::notdone() { return ( i < dimd ); }


////////////////////////////////////////////////////////////////
// a pointer that loops over all elements
// of an Idx

// empty constructor;
template <class T> IdxIter<T>::IdxIter() { };

template <class T> T *IdxIter<T>::init(Idx<T> &idx) {
  iterand = &idx;
  i = 0;
  j = iterand->spec.ndim;
  data = iterand->storage->data + iterand->spec.offset;
  if (iterand->spec.contiguousp()) {
    d[0] = -1;
  } else {
    for(int i=0; i < iterand->spec.ndim; i++) { d[i] = 0; }
  }
  n = iterand->spec.nelements();
  return data;
}

template <class T> T *IdxIter<T>::next() {
  i++;
  if (d[0] < 0) {
    // contiguous Idx
    data++;
  } else {
    // non-contiguous Idx
    j--;
    do {
      if (j<0) break;
      if (++d[j] < iterand->spec.dim[j]) {
	data += iterand->spec.mod[j];
	j++;
      } else {
	data -= iterand->spec.dim[j] * iterand->spec.mod[j];
	d[j--] = -1;
      }
    } while (j < iterand->spec.ndim);
  }
  return data;
}

template <class T> bool IdxIter<T>::notdone() { return ( i < n ); }

#endif // IF USING_STL_ITERS == 1

} // namespace ebl

#endif