idx.cpp

Gaussian Mixture Algorithm

    if (s5>=0) { 
      if (!ndimset) { setndim(6); ndimset = true; } 
      dim[5] = s5; mod[5] = md; md *= s5; 
    } else { if (ndimset) { throw(-6); } }
    if (s4>=0) { 
      if (!ndimset) { setndim(5); ndimset = true; } 
      dim[4] = s4; mod[4] = md; md *= s4; 
    } else { if (ndimset) { throw(-5); } }
    if (s3>=0) { 
      if (!ndimset) { setndim(4); ndimset = true; } 
      dim[3] = s3; mod[3] = md; md *= s3; 
    } else { if (ndimset) { throw(-4); } }
    if (s2>=0) { 
      if (!ndimset) { setndim(3); ndimset = true; } 
      dim[2] = s2; mod[2] = md; md *= s2; 
    } else { if (ndimset) { throw(-3); } }
    if (s1>=0) { 
      if (!ndimset) { setndim(2); ndimset = true; } 
      dim[1] = s1; mod[1] = md; md *= s1; 
    } else { if (ndimset) { throw(-2); } }
    if (s0>=0) { 
      if (!ndimset) { setndim(1); ndimset = true; } 
      dim[0] = s0; mod[0] = md; md *= s0; 
    } else { if (ndimset) { throw(-1); } }
    if (!ndimset) { setndim(0); }
  }
  catch(int v) { ylerror("IdxSpec: bad dimensions in constructor"); }
}

// generic constructor for any dimension.
// The dim and mod arrays past as argument are copied.
IdxSpec::IdxSpec(intg o, int n, intg *ldim, intg *lmod) {
  DEBUG("IdxSpec::IdxSpec: %ld\n",(intg)this);
  dim = NULL; mod = NULL;
  offset = o;
  setndim(n);
  for (int i=0; i<n; i++) { 
    if ( dim[i]<0 ) { ylerror("negative dimension"); }
    dim[i] = ldim[i]; mod[i] = lmod[i]; }
}



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

// pretty print 
void IdxSpec::pretty(FILE *f) {
  int i;
  fprintf(f,"  IdxSpec %ld\n",(intg)this);
  fprintf(f,"    ndim=%d\n",ndim);
  fprintf(f,"    offset=%ld\n",offset);
  if (ndim>0) {
    fprintf(f,"    dim=[ "); 
    for (i=0; i<ndim-1; i++){ fprintf(f,"%ld, ",dim[i]); }
    fprintf(f,"%ld]\n",dim[ndim-1]); 
    fprintf(f,"    mod=[ "); 
    for (i=0; i<ndim-1; i++){ fprintf(f,"%ld, ",mod[i]); }
    fprintf(f,"%ld]\n",mod[ndim-1]);
  } else {
    fprintf(f,"    dim = %ld, mod = %ld\n",(intg)dim, (intg)mod);
  }
  fprintf(f,"    footprint= %ld\n",footprint());
  fprintf(f,"    contiguous= %s\n",(contiguousp())?"yes":"no");
}

void IdxSpec::pretty(std::ostream& out) {
  int i;
  out << "  IdxSpec " << (intg)this << "\n";
  out << "    ndim= " << ndim << "\n";
  out << "    offset= " << offset << "\n";
  if (ndim>0) {
    out << "    dim=[ ";
    for (i=0; i<ndim-1; i++){ out << dim[i] << ", "; }
    out << dim[ndim-1] << "]\n"; 
    out << "    mod=[ "; 
    for (i=0; i<ndim-1; i++){ out << mod[i] << ", "; }
    out << mod[ndim-1] << "]\n";
  } else {
    out << "    dim = " << (intg)dim << ", mod = " << (intg)mod <<"\n";
  }
  out << "    footprint= " << footprint() << "\n";
  out << "    contiguous= " << ((contiguousp())? "yes":"no") << "\n";
}
////////////////////////////////////////////////////////////////
// select, narrow, unfold, etc
// Each function has 3 version: 
// 1. XXX_into: which writes the result
// into an existing IdxSpec apssed as argument.
// 2. XXX_inplace: writes into the current IdxSpec
// 3. XXX: creates a new IdxSpec and returns it.

intg IdxSpec::select_into(IdxSpec *dst, int d, intg n) {
  try {
    if (ndim <= 0) throw("cannot select a scalar");
  }
  catch(const char *s) { ylerror(s); return -1;}
  // this preserves the dim/mod arrays if dst == this
  dst->setndim(ndim-1);
  dst->offset = offset + n * mod[d];
  if (ndim -1 > 0) { // dim and mod don't exist for Idx0
  	for (int j=0; j<d; j++) {
  		dst->dim[j] = dim[j];
  		dst->mod[j] = mod[j];
  	}
  	for (int j=d; j<ndim-1; j++) {
  		dst->dim[j] = dim[j+1];
  		dst->mod[j] = mod[j+1];
  	}
  }
  return n;
}

intg IdxSpec::select_inplace(int d, intg n) {
  return select_into(this, d, n);
}

IdxSpec IdxSpec::select(int d, intg n) {
  // create new IdxSpec of order ndim-1
  IdxSpec r;
  select_into(&r, d, n);
  return r;
}
  
////////////////////////////////////////////////////////////////

intg IdxSpec::narrow_into(IdxSpec *dst, int d, intg s, intg o) {
  try {
    if (ndim <= 0) throw("cannot narrow a scalar");
    if ((d < 0) || (d>=ndim)) throw("narrow: illegal dimension index");
    if ((o < 0)||(s < 1)||(s+o > dim[d])) throw("narrow: illegal size/offset");
  }
  catch(const char *s) { ylerror(s); return -1;}
  // this preserves the dim/mod arrays if dst == this
  dst->setndim(ndim);
  dst->offset = offset + o * mod[d];
  for (int j=0; j<ndim; j++) {
    dst->dim[j] = dim[j];
    dst->mod[j] = mod[j];
  }
  dst->dim[d] = s;
  return s;
}

intg IdxSpec::narrow_inplace(int d, intg s, intg o) {
  return narrow_into(this, d, s, o);
}

IdxSpec IdxSpec::narrow(int d, intg s, intg o) {
  // create new IdxSpec of order ndim
  IdxSpec r;
  narrow_into(&r, d, s, o);
  return r;
}

////////////////////////////////////////////////////////////////
// transpose

// tranpose two dimensions into pre-existing IdxSpec
int IdxSpec::transpose_into(IdxSpec *dst, int d1, int d2) {
  try {
    if ((d1 < 0) || (d1 >= ndim) || 
	(d2 < 0) || (d2 >= ndim)) 
      throw("tranpose: illegal dimension index");
  }
  catch(const char *s) { ylerror(s); return -1;}
  // this preserves the dim/mod arrays if dst == this
  dst->setndim(ndim);
  dst->offset = offset;
  for (int j=0; j<ndim; j++) {
    dst->dim[j] = dim[j];
    dst->mod[j] = mod[j];
  }
  intg tmp;
  // we do this in case dst = this
  tmp=dim[d1]; dst->dim[d1]=dim[d2]; dst->dim[d2]=tmp;
  tmp=mod[d1]; dst->mod[d1]=mod[d2]; dst->mod[d2]=tmp;
  return ndim;
}

// tranpose all dims with a permutation vector
int IdxSpec::transpose_into(IdxSpec *dst, int *p) {
  try {
    for (int i=0; i<ndim; i++) {
      if ((p[i] < 0) || (p[i] >= ndim)) throw("tranpose: illegal dimension index");
    }
  }
  catch(const char *s) { ylerror(s); return -1;}
  dst->setndim(ndim);
  dst->offset = offset;
  if (dst == this) {
    // we need temp storage if done in place
    intg tmpdim[MAXDIMS], tmpmod[MAXDIMS];
    for (int j=0; j<ndim; j++) {
      tmpdim[j] = dim[p[j]];
      tmpmod[j] = mod[p[j]];
    }
    for (int j=0; j<ndim; j++) {
      dst->dim[j] = tmpdim[j];
      dst->mod[j] = tmpdim[j];
    }
  } else {
    // not in place
    for (int j=0; j<ndim; j++) {
      dst->dim[j] = dim[p[j]];
      dst->mod[j] = mod[p[j]];
    }
  }
  return ndim;
}


int IdxSpec::transpose_inplace(int d1, int d2) {
  return transpose_into(this, d1, d2);
}

int IdxSpec::transpose_inplace(int *p) {
  return transpose_into(this, p);
}

IdxSpec IdxSpec::transpose(int d1, int d2) {
  IdxSpec r;
  transpose_into(&r, d1, d2);
  return r;
}

IdxSpec IdxSpec::transpose(int *p) {
  IdxSpec r;
  transpose_into(&r, p);
  return r;
}


////////////////////////////////////////////////////////////////
// unfold

// d: dimension; k: kernel size; s: stride.
intg IdxSpec::unfold_into(IdxSpec *dst, int d, intg k, intg s) {
  intg ns; // size of newly created dimension
  try {
    if (ndim <= 0) throw("cannot unfold an Idx of maximum order");
    if ((d < 0) || (d>=ndim)) throw("unfold: illegal dimension index");
    if ((k < 1) || (s < 1)) throw("unfold: kernel and stride must be >= 1");
    ns = 1+ (dim[d]-k)/s;
    if ((ns <= 0) || ( dim[d] != s*(ns-1)+k )) throw("unfold: kernel and stride incompatible with size");
  }
  catch(const char *s) { ylerror(s); return -1;}
  // this preserves the dim/mod arrays if dst == this
  dst->setndim(ndim+1);
  dst->offset = offset;
  for (int i=0; i<ndim; i++) {
    dst->dim[i] = dim[i];
    dst->mod[i] = mod[i];
  }
  dst->dim[ndim] = k;
  dst->mod[ndim] = mod[d];
  dst->dim[d] = ns;
  dst->mod[d] = mod[d]*s;
  return ns;
}

intg IdxSpec::unfold_inplace(int d, intg k, intg s) {
  return unfold_into(this, d, k, s);
}

IdxSpec IdxSpec::unfold(int d, intg k, intg s) {
  IdxSpec r;
  unfold_into(&r, d, k, s);
  return r;
}

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

// return true if two idxspec have the same dimensions,
// i.e. if all their dimensions are equal (regardless of strides).
bool same_dim(IdxSpec &s1, IdxSpec &s2) {
  if ( s1.ndim != s2.ndim ) return false; 
  for (int i=0; i<s1.ndim; i++) { if (s1.dim[i] != s2.dim[i]) return false; }
  return true;
}

} // end namespace ebl