/**
Allocate variable vector of matrices with dimensions
[sl to sh] x ([nrl to nrh] x [ncl to nch])
where nrl is a vector of indexes.
\param sl lower index of vector
\param sh upper index of vector
\param nrl vector of lower row indexes for matrix
\param nrh upper row index for matrix
\param ncl upper column index for matrix
\param nch upper column index for matrix
*/
void dvar3_array::allocate(int sl, int sh, const ivector& nrl, int nrh,
  int ncl, int nch)
{
  if (sh < sl)
  {
    allocate();
    return;
  }
  if (sl !=nrl.indexmin() || sh !=nrl.indexmax())
  {
    cerr << "Incompatible array bounds in "
     "dmatrix(int nrl,int nrh, const ivector& ncl, const ivector& nch)" << endl;
    ad_exit(1);
  }
  if ( (shape=new three_array_shape(sl,sh)) == 0)
  {
    cerr << " Error allocating memory in dvar3_array contructor" << endl;
  }
  if ( (t = new dvar_matrix[slicesize()]) == 0)
  {
    cerr << " Error allocating memory in dvar3_array contructor" << endl;
    ad_exit(21);
  }
  t -= slicemin();
  for (int i = sl; i <= sh; ++i)
  {
    t[i].allocate(nrl(i), nrh, ncl, nch);
  }
}

  bool MLP::calcGradient(const dmatrix& inputs,
                         const ivector& ids,
                         dvector& grad) {

    if (inputs.rows() != ids.size()) {
      setStatusString("Number of vectors not consistent with number of ids");
      return false;
    }

    dvector tmp;
    int i;
    double tmpError;

    totalError = 0;
    calcGradient(inputs.getRow(0),ids.at(0),grad);
    computeActualError(ids.at(0),totalError);

    for (i=1;i<inputs.rows();++i) {
      calcGradient(inputs.getRow(i),ids.at(i),tmp);
      computeActualError(ids.at(i),tmpError);
      grad.add(tmp);
      totalError+=tmpError;
    }

    return true;
  }

  //get most(best) available label from histogram
  int kNearestNeighFilter::getMostLabel(const ivector& histogram,
			       const imatrix& src,
			       const int& row, const int& col) const{

    int numOfMax = 0;
    int maxIndex = -1; // first index, which is max
    int max = 0; //
    for(int i=0;i<histoSize;++i) {
      if(histogram.at(i) < max); // for speed up (probability)
      else if(histogram.at(i) > max) {
	max = histogram.at(i);
	numOfMax = 1;
	maxIndex = i;
      }
      else //if(histogram.at(i) == max)
	++numOfMax;
    }

    //is there more than one possibility ?
    if (numOfMax == 1)
      return maxIndex;
    // is the kernel center one of the max's?
    else if(histogram.at(src.at(row,col)) == max)
      return src.at(row,col);
    else
      return getMedian(histogram,max,numOfMax);
  };

/**
 * Description not yet available.
 * \param
 */
 void dmatrix::allocate(int nrl,int nrh,const ivector& ncl,int nch)
 {
   if (nrh<nrl)
   {
     allocate();
     return;
   }
   if (nrl !=ncl.indexmin() || nrh !=ncl.indexmax())
   {
     cerr << "Incompatible array bounds in dmatrix(int nrl,int nrh,"
       " int ncl,const ivector& nch)\n";
     ad_exit(1);
   }
   index_min=nrl;
   index_max=nrh;
   if ( (m = new dvector [rowsize()]) == 0)
   {
     cerr << " Error allocating memory in dmatrix contructor\n";
     ad_exit(21);
   }
   if ( (shape = new mat_shapex(m))== 0)
   {
     cerr << " Error allocating memory in dmatrix contructor\n";
     ad_exit(21);
   }
   m -= rowmin();
   for (int i=rowmin(); i<=rowmax(); i++)
   {
     m[i].allocate(ncl(i),nch);
   }
 }

/**
 * Description not yet available.
 * \param
 */
 void lmatrix::allocate(int nrl, int nrh, int ncl, const ivector& nch)
 {
   if (nrl !=nch.indexmin() || nrh !=nch.indexmax())
   {
     cerr << "Incompatible array bounds in "
     "lmatrix::allocate(int nrl,int nrh,int ncl, const ivector& nch)\n";
     ad_exit(1);
   }
   if ( (shape = new mat_shape(nrl,nrh,ncl,nch(nch.indexmin())))== 0)
   {
     cerr << " Error allocating memory in lmatrix contructor\n";
     ad_exit(21);
   }
   size_t rs=rowsize();
   if ( (m = new lvector [rs]) == 0)
   {
     cerr << " Error allocating memory in lmatrix contructor\n";
     ad_exit(21);
   }
   m -= rowmin();
   for (int i=rowmin(); i<=rowmax(); i++)
   {
     m[i].allocate(ncl,nch(i));
   }
 }

/**
 * Description not yet available.
 * \param
 */
dvector::dvector(const ivector& u)
 {
   allocate(u.indexmin(),u.indexmax());
   for ( int i=indexmin(); i<=indexmax(); i++)
   {
     elem(i)=u.elem(i);
   }
 }

/**
Returns the sum of the squares of all elements in ivec.

\param ivec ivector
*/
int norm2(const ivector& ivec)
{
  int sum = 0;
  for (int i = ivec.indexmin(); i <= ivec.indexmax(); ++i)
  {
    sum += ivec(i) * ivec(i);
  }
  return sum;
}

/**
 * Description not yet available.
 * \param
 */
dvector dvector::operator ()(const ivector& u)
 {
   dvector tmp(u.indexmin(),u.indexmax());

   for ( int i=u.indexmin(); i<=u.indexmax(); i++)
   {
     tmp(i)=(*this)(u(i));
   }
   return tmp;
 }

/**
Return dvector results of squaring elements in a values;
constant vector object.

\ingroup misc
\param values of constant object to be squared.
\return vector of the same length as #values containing \f$values_i^2\f$
*/
ivector square(const ivector& values)
{
  ivector results;
  results.allocate(values);
  for (int i = values.indexmin(); i <= values.indexmax(); ++i)
  {
    results(i) = values(i) * values(i);
  }
  return results;
}

  // the kernel runs inside the image
  void kNearestNeighFilter::histogramMethodMiddle(const imatrix& src,
                                              imatrix& dest,
                                              ivector& histogram,
                                              const int& row,int& col) const {
    int i,j;//index
    int numOfMax, maxIndex;
    int max=0;
    const int maxChange = sizeOfKernel+1;//max change for "max"

    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int lastCol = src.lastColumn()-limit;
    const int r = row+limit;
    col = limit;

    int v; //del test
    while(col <= (lastCol-1)) {
      j = col-limit;
      // sub labels left form the kernel
      for(i=row-limit;i<=r;++i) {
        --histogram.at(src.at(i,j));
      }
      // add labels right from the kernel
      ++col;
      j = col+limit;
      for(i=row-limit;i<=r;++i) {
	v = src.at(i,j);
        ++histogram.at(src.at(i,j));
      }

      //get most(best) available label
      numOfMax = 0;
      maxIndex = -1;
      max -= maxChange; //=0;
      for(i=0;i<histoSize;++i) {
	if(histogram.at(i) < max);// for speed up (probability)
	else if(histogram.at(i) > max) {
	  max = histogram.at(i);
	  numOfMax = 1;
	  maxIndex = i;
	}
	else //if(histogram.at(i) == max)
	  ++numOfMax;
      }

      //is there more than one possibility ?
      if(numOfMax == 1)
	dest.at(row,col) =  maxIndex;
      // is the kernel center one of the max's?
      else if(histogram.at(src.at(row,col)) == max)
	dest.at(row,col) = src.at(row,col);
      else
	dest.at(row,col) = getMedian(histogram,max,numOfMax);
    }//while
  };

 void allocate(int lb, int ub)
 {
     indx.allocate(lb, ub);
     indx2.allocate(lb, ub);
     ivector iv(lb + 1, ub);
     iv.fill_seqadd(lb, 1);
     L.allocate(lb + 1, ub, lb, iv);
     ivector iv1(lb, ub);
     iv1.fill_seqadd(lb, 1);
     U.allocate(lb, ub, lb, iv1);
     indx2.fill_seqadd(lb, 1);
 }

 void svm::makeTargets(const ivector& ids) {
   // expand each class label i to a vector v with v[j]=1 if j == i,
   // and j[j]=-1 if j != i
   srcIds=ids;
   dmatrix* t=new dmatrix(nClasses,ids.size(),-1.0);
   // iterate over training labels
   for (int i=0; i<t->columns(); i++) {
     t->at(idMap[ids.at(i)],i)=1;
   }
   if (target != 0) {
     delete target;
   }
   target=t;
 }

  int kNearestNeighFilter::getMedian(const ivector& histogram,
                                     const int max,
                                     const int numOfMax) const {
    ivector vect(numOfMax,0);
    int i,z=0;
    const int size=histogram.size();
    for(i=0;i<size;++i) {
      if (histogram.at(i) == max) {
        vect.at(z++) = i;
      }
    }

    return vect.at(z/2);
  }

  void svm::buildIdMaps(const ivector& ids) {
    int j=0;
    // create reverse id map
    idMap.clear();
    for (int i=0; i<ids.size(); i++) {
      if (idMap.find(ids.at(i)) == idMap.end()) {
      _lti_debug("Mapping external id " << ids.at(i) << " to " << j << std::endl);
        rIdMap[j]=ids.at(i);
        idMap[ids.at(i)]=j++;
      }
    }

    nClasses=j;
  }

  void MLP::checkHowManyOutputs(const ivector& ids) {
    // count how many different ids are present in the training set
    std::map<int,int> extToInt;
    std::map<int,int>::iterator it;

    int i,k;
    for (i=0,k=0;i<ids.size();++i) {
      it = extToInt.find(ids.at(i));
      if (it == extToInt.end()) {
        extToInt[ids.at(i)] = k;
        ++k;
      }
    }

    outputs = extToInt.size();
  }

void CKDTree::SearchTree(KD_Node *parent, Point &p, float radius, ivector &vec_r)
{
	if (parent==NULL || parent->IsSearched)
		return;
	parent->IsSearched=true;
	//is a child node
	if(parent->m_left==NULL && parent->m_right==NULL)
	{
		int nSize=parent->m_vecIn.size();
		for (int i=0; i<nSize; ++i)
			if(m_regionArray[parent->m_vecIn[i]].IsSeen(p, radius))
				vec_r.push_back(parent->m_vecIn[i]);
	}
	//not a child node
	else
	{
		if(parent->m_left!=NULL && !parent->m_left->IsSearched && 
			(parent->m_left->m_regOverlap.IsSeen(p, radius)))
			SearchTree(parent->m_left, p, radius, vec_r);
		if (parent->m_right!=NULL && !parent->m_right->IsSearched &&
			(parent->m_right->m_regOverlap.IsSeen(p, radius)))
			SearchTree(parent->m_right, p, radius, vec_r);
	}
	//back trace : search parent
	SearchTree(parent->m_parent, p, radius, vec_r);
}

/**
Set elements of ivec to zero starting from level + 1;

\param level is the index of ivec
*/
void clean(ivector& v, int level)
{
  int max = v.indexmax();
  for (int i = level + 1; i <= max; ++i)
  {
    v(i) = 0;
  }
}

 void pvm_int::assign(const ivector& u)
 {
   if(ad_comm::pvm_manager)
   {
     int nsp=ad_comm::pvm_manager->num_slave_processes;
     if (u.indexmin() !=0 || u.indexmax() != nsp)
     {
       cerr << "Error in pvm_int::assign  validindex bounds must be 0 "
            << ad_comm::pvm_manager->num_slave_processes << endl;
       ad_exit(1);
     }
     if (allocated(v))
       v.deallocate();
     v.allocate(0,nsp);
     v=u;
     d=u(0);
   }
 }

 cltudecomp(int lb, int ub):indx(lb, ub), indx2(lb, ub)
 {
     ivector iv(lb + 1, ub);
     iv.fill_seqadd(lb, 1);
     L.allocate(lb + 1, ub, lb, iv);
     ivector iv1(lb, ub);
     iv1.fill_seqadd(lb, 1);
     U.allocate(lb, ub, lb, iv1);
     indx2.fill_seqadd(lb, 1);
 }

  // return probability value of an rgb pixel
  float probabilityMap2D::apply(const ubyte &value1, const ubyte &value2, ivector& theBin) const {

    assert((probabilityHistogram.dimensions() == 2) &&
           (theBin.size() == 2));
    
    theBin[0] = lookupTable[0][static_cast<int>(value1)];
    theBin[1] = lookupTable[1][static_cast<int>(value2)];

    return static_cast<float>(probabilityHistogram.at(theBin));
  }