#include "multiDS.h"


#include "../Matrices/binMatrix.h"
#include "../Matrices/setCovering.h"
#include "../LP/linearProgram.h"
#include "main_def.h"
								   


#define _CHECK_INDEX



//////////////////
// Constructors //
//////////////////

tcT multiDS<T>::multiDS () : 
   ds(0), label(), nbSets(0), list(), normalized(false) {}

tcT multiDS<T>::multiDS (int nbCategories) :
   ds(new dataSet<T> [nbCategories]), label(nbCategories), 
   nbSets(nbCategories), list(), normalized(false) {}


tcT multiDS<T>::multiDS (const multiDS<T>& from, float rate) :
   ds(new dataSet<T> [from.nbCats()]), label(from.label),
   nbSets(from.nbSets), list(), normalized(from.normalized)
{
   for (int set=0; set<nbCats(); set++) 
      ds[set] = dataSet<T>(from.ds[set], rate);
}


tcT multiDS<T>::multiDS (const multiDS<T>& from, const Array<int>& varIndex) :
   ds(new dataSet<T> [from.nbCats()]), label(from.label),
   nbSets(from.nbSets), list(), normalized(false)
{
   for (int set=0; set<nbCats(); set++)
      ds[set] = dataSet<T>(from.ds[set], varIndex);
}


tcT multiDS<T>::multiDS (const dataSet<T>& from) :
   ds(0), label(), nbSets(0), list(), normalized(false)
{	
#ifdef _CHECK_INDEX
   if ( from.classCol()<0 || from.classCol()>=from.dim() )
      throw whereOutOfRange(
	 "multiDS<T>::multiDS(const dataSet<T>&, int from.classCol())",
	 "from.classCol()",from.classCol(),0,from.dim()-1);
#endif
   // determine the number of distinct categories
   Matrix<int> singleton(1,1,from.classCol()); int garbage;
   dataSet<T> categories(from,singleton);
   categories.sort().checkMult(garbage);
   Matrix<T> catIndex(categories.matrix().ta());
   ds = new dataSet<T> [catIndex.n()]; label.resize(catIndex.n());
   for (int cI=0; cI<catIndex.n(); cI++)
   {
      Matrix<T> cm(from.distinct(),from.dim()-1);
      Matrix<int> multiple(from.distinct());
      Matrix<int> labelle(from.distinct());
      int cmPtr=0;
      for (int dPtr=0; dPtr<from.distinct(); dPtr++)
	 if ( from(dPtr,from.classCol())==catIndex(cI) )
	 {
	    if ( from.classCol() > 0 )
	       cm.s(cmPtr,1,0,from.classCol()) 
	       = from(dPtr).s(0,1,0,from.classCol());
	    if ( from.classCol() < from.dim()-1 )
	       cm.s(cmPtr,1,from.classCol()) 
	       = from(dPtr).s(0,1,from.classCol()+1);
	    multiple(cmPtr) = from.multiplicity(dPtr);
	    labelle(cmPtr) = from.label(dPtr);
	    cmPtr++;
	 }
      if ( cmPtr )
      {
	 ds[nbSets] = dataSet<T>(	cm.resize_(cmPtr,cm.n()),
				multiple.resize_(cmPtr),
				labelle.resize_(cmPtr));
	 label(nbSets) = catIndex(cI);
	 nbSets++;
      }
   }
   if ( nbSets < catIndex.n() ) // some categories were inexistant
   {
      label.resize_(nbSets);
      dataSet<T>* aux = ds;
      ds = new dataSet<T> [nbSets];
      for (int i=0; i<nbSets; i++) ds[i] = aux[i];
      delete [] aux;
   }
   for (int set=0; set<nbSets; set++)
      ds[set].monotone.resize(from.monotone);
}


tcT multiDS<T>::multiDS (	const dataSet<T>& from, const Matrix<int>& catIndex) :
   ds(new dataSet<T> [catIndex.n()]), label(catIndex.n()), 
   nbSets(0), list(), normalized(false)
{
#ifdef _CHECK_INDEX
   if ( from.classCol()<0 || from.classCol()>=from.dim() )
      throw whereOutOfRange(
	 "multiDS<T>::multiDS(const dataSet<T>&, int, const Matrix<T>&)",
	 "from.classCol()",from.classCol(),0,from.dim());
#endif
   for (int cI=0; cI<catIndex.n(); cI++)
   {
      Matrix<T> cm(from.distinct(),from.dim()-1);
      Matrix<int> multiple(from.distinct());
      Matrix<int> labelle(from.distinct());
      int cmPtr=0;
      for (int dPtr=0; dPtr<from.distinct(); dPtr++)
	 if ( from(dPtr,from.classCol())==catIndex(cI) )
	 {
	    if ( from.classCol() > 0 )
	       cm.s(cmPtr,1,0,from.classCol()) 
	       = from(dPtr).s(0,1,0,from.classCol());
	    if ( from.classCol() < from.dim()-1 )
	       cm.s(cmPtr,1,from.classCol()) 
	       = from(dPtr).s(0,1,from.classCol()+1);
	    multiple(cmPtr) = from.multiplicity(dPtr);
	    labelle(cmPtr) = from.label(dPtr);
	    cmPtr++;
	 }
      if ( cmPtr )
      {
	 ds[nbSets] = dataSet<T>(	cm.resize_(cmPtr,cm.n()),
				multiple.resize_(cmPtr),
				labelle.resize_(cmPtr));
	 label(nbSets) = catIndex(cI);
	 nbSets++;
      }
   }
   if ( nbSets++ < catIndex.n() ) // some categories were inexistant
   {
      label.resize_(nbSets);
      dataSet<T>* aux = ds;
      ds = new dataSet<T> [nbSets];
      for (int i=0; i<nbSets; i++) ds[i] = aux[i];
      delete [] aux;
   }
   for (int set=0; set<nbSets; set++)
      ds[set].monotone.resize(from.monotone);
}


tcT multiDS<T>::~multiDS ()
{
   delete [] ds;
}



///////////////
// Operators //
///////////////

tcT /*inline*/ int multiDS<T>::nbCats() const { return nbSets; }

tcT /*inline*/ int multiDS<T>::dim() const { return ds[0].dim(); }


tcT int multiDS<T>::nbPairs() const 
{
   int nbConstrs=0;
   for (int clL=0; clL<nbCats()-1; clL++)
      for (int clR=clL+1; clR<nbCats(); clR++)
	 nbConstrs += (ds[clL].distinct() * ds[clR].distinct());
   return nbConstrs;
}


tcT int multiDS<T>::distinct() const
{
   int res=0;
   for (int cl=0; cl<nbCats(); cl++) res += ds[cl].distinct();
   return res;
}		


tcT int multiDS<T>::card() const
{
   int res=0;
   for (int cl=0; cl<nbCats(); cl++) res += ds[cl].card();
   return res;
}


tcT /*inline*/ char multiDS<T>::monotone(const int attr) const
{
   if ( !ds[0].monotone.n() )
      return nonMonotonicAttr;
   else
      return ds[0].monotone(attr);
}


tcT multiDS<T>& multiDS<T>::checkMult()
{
   int saved; 
   for (int set=0; set<nbCats(); set++) ds[set].checkMult(saved);
   return *this;
}


tcT multiDS<T>& multiDS<T>::sort()
{
   for (int set=0; set<nbCats(); set++) ds[set].sort();
   return *this;
}


tcT multiDS<T>& multiDS<T>::reorganize()
{
   for (int set=0; set<nbCats(); set++) ds[set].reorganize(true);
   return *this;
}


template <class T>
ostream& operator << (ostream& s, const multiDS<T>& myself)
{
   for (int set=0; set<myself.nbCats(); set++)
   {
      s  << "\n----------c-l-a-s-s--- " << myself.label(set)
	 << " ----------\n" << myself.ds[set];
   }
   return s << endl;
}


tcT void multiDS<T>::split(multiDS<T>& first, multiDS<T>& second, float rate) const
{
   if ( first.nbCats()!=nbCats() ) 
   {
      delete first.ds; first.ds = new dataSet<T> [nbCats()];
      first.nbSets=nbCats();
      first.list.resize();
   }
   if ( second.nbCats()!=nbCats()) 
   {
      delete second.ds; second.ds = new dataSet<T> [nbCats()];
      second.nbSets=nbCats();
      second.list.resize();
   }
   first.label.resize(label); first.normalized=normalized;
   second.label.resize(label); second.normalized=normalized;
   for (int set=0; set<nbCats(); set++)
   {
      ds[set].split(first.ds[set],second.ds[set],rate);
   }
   first.list.resize(0); second.list.resize(0);
}

tcT void multiDS<T>::partition(Matrix<int>& part, int K) const
   // part is a maxCard+K+1 x nbCats() matrix. 
   // part(0,s)=K for all s.
   // part(k,s)=cardinality of fold k in set s, k=1,...,K.
   // part(K+i,s)=1,...,K, label of the i^th el. of set s, i=1,...,ds[s].card().
{
   // compute the maximum card among all datasets
   int maxCard=0;
   for (int set=0; set<nbCats(); set++)
      maximize(maxCard,ds[set].card());
   part.resize(maxCard+K+1,nbCats());
   part=0; part.s(0)=K;

   // compute the cardinalities of each set in each part.
   // In a first phase, part(k,s) will be set to floor(ds[s].card()/K), 
   // and then, for some of the k, cardinal(k,s) will be increased by at 
   // most 1, until sum_k part(k,s) = ds[s].card().
   // Moreover, this will be done in such a way that there will be never two
   // parts k1 and k2 such that sum_s part(k1,s) and sum_s part(k2,s)
   // differ from more than 1.

   Matrix<int> smallerPart(1,K,1); 
   int nbSmaller = K;
   // For the purpose of balancing all the parts, the number of elements
   // in two different parts will never differ from more than 1. smallerPart(k)
   // indicates those parts with a currently smaller number of elements and
   // nbSmaller = smallerPart.sum().

   for (set=0; set<nbCats(); set++)
   {
      // cardinalities of all parts into set set is set to the floor values.
      const int theFloor = floor(float(ds[set].card())/K);
      part.s(1,K,set,1) = theFloor;

      int toAdd = ds[set].card() - K * theFloor;
      if ( toAdd>0 && nbSmaller<=toAdd )
      {
	 part.s(1,K,set,1) += smallerPart.t();
	 toAdd -= nbSmaller;
	 nbSmaller = K; smallerPart = 1;
      }
      if ( toAdd>0 )
      {
	 // choose at random toAdd parts among the nbSmaller ones for which
	 // smallerPart is 1.
	 int skipped=0;
	 for (int k=0; k<K; k++)
	    if ( smallerPart(k) && part(k+1,set)==theFloor &&
		 randomR() < float(toAdd)/(nbSmaller - skipped++) )
	    {
	       part(k+1,set)++;
	       smallerPart(k)=0; nbSmaller--; skipped--;
	       if ( !(--toAdd) ) break;
	    }
      }

      // update PART such that elements of SET are shuffled through all the
      // parts according to the cardinalities.
      int finger=0;
      for (int p=1; p<=K; finger+=part(p,set), p++)
	 part.s(K+1+finger,part(p,set),set,1) = p;

      // reorder at random the partition
      for (int data=0; data<ds[set].card(); data++)
	 Exchange(part(K+1+data,set),part(K+1+randomI(ds[set].card()),set));
   }
}


tcT void multiDS<T>::setFold(multiDS<T>& data, 
			   const Matrix<int>& part, const int fold) const
{
   if ( data.nbCats()!=nbCats() ) 
   {
      delete data.ds; data.ds = new dataSet<T> [nbCats()];
      data.nbSets=nbCats();
      data.list.resize();
   }
   const int K = part(0,0);
   data.label.resize(label); data.normalized=normalized;
   for (int set=0; set<nbCats(); set++)
   {
      const int cardin = 
	 (fold>0 ? part(fold,set) : ds[set].card()-part(-fold,set));
      Matrix<T> d(cardin,dim());
      Matrix<int> multiple(1,cardin,1);
      Matrix<int> labelle(1,cardin);
      int ptrD=0; int ptrDS=0; int mu=1;
      for (int ptrP=0; ptrP<ds[set].card(); ptrP++, mu++)
      {
	 if ( mu > ds[set].multiplicity(ptrDS) )
	 {
	    mu=1;
	    ptrDS++;
	 }
	 if ( fold>0 ? part(K+1+ptrP,set)==fold : part(K+1+ptrP,set)!=(-fold) )
	 {
	    d.s(ptrD)=ds[set](ptrDS);
	    labelle(ptrD)=ds[set].label(ptrDS);
	    ptrD++;
	 }
      }
      data.ds[set] = dataSet<T>(d,multiple,labelle,ds[set].monotone);
   }
}


tcT boolean multiDS<T>::isConsistant(const boolean signalIfNot/**=true**/) const
{
   boolean res=true;
   for (int set1=0; set1<nbCats()-1; set1++)
      for (int set2=set1+1; set2<nbCats(); set2++)
	 for (int data1=0; data1<ds[set1].distinct(); data1++)
	    for (int data2=0; data2<ds[set2].distinct(); data2++)
	    {
	       boolean different=false;
	       for (int attr=0; attr<dim(); attr++)
		  if ( ( ds[set1](data1,attr)<ds[set2](data2,attr) &&
		         monotone(attr)!=negativeAttr ) ||
		       ( ds[set1](data1,attr)>ds[set2](data2,attr) &&
		         monotone(attr)!=positiveAttr )
		     )
		  { different=true; break; }
	       if ( !different )
		  if ( signalIfNot )
		  {
		     res=false;
		     flog<<"\n### Data" << setw(6) << ds[set1].label(data1)
		         << " cannot be distinguished from data" << setw(6) 
			 << ds[set2].label(data2) << endl
			 << ds[set1](data1) << ds[set2](data2) ;
		  }
		  else return false;
	    }
   return res;
}


tcT boolean isInBox
(Array<T>& corner1, Array<T>& corner2, Array<T>& candidate)
{
//	Check whether candidate[0] is in the box delimitated by the two
   for (int v = 1; v < corner1.n(); v++)
      // v start at one since 0 is already OK due to the sorting
      if ( corner1(v) < corner2(v) )
      {
	 if ( candidate(v) < corner1(v) ||
	      candidate(v) > corner2(v) )
	    return false;
      }
      else if ( candidate(v) < corner2(v) ||
		candidate(v) > corner1(v) )
	 return false;
   return true;
}


tcT Array<int> multiDS<T>::listOfPairs(
   boolean compress/**=false**/, boolean recreate/**=false**/)
{
   if ( !recreate && list.n() ) return list;
   int clL, clR;
   list.resize(nbPairs(),4); 
   int ptr=0; 
   boolean addThis=true;
   for (clL=0; clL<nbCats()-1; clL++)
      for (clR=clL+1; clR<nbCats(); clR++)
	 for (int l=0; l<ds[clL].distinct(); l++)
	 {
	    boolean l_complete=true;
	    int i;
	    for(i=0;i< dim();i++)
	       if (unknown(ds[clL](l,i))) {l_complete=false; break;}
	    for (int r=0; r<ds[clR].distinct(); r++)
	    {
	       boolean r_complete=true;
	       if (l_complete)
		  for(i=0;i< dim();i++)
		     if (unknown(ds[clR](r,i))) {r_complete=false; break;}
	       addThis=true;
	       if ( compress && l_complete && r_complete) 
		  // applies the box's idea for suppressing some pairs
	       {
		  T min0 = minimum(ds[clL](l,0),ds[clR](r,0));
		  T max0 = maximum(ds[clL](l,0),ds[clR](r,0));
		  for (int clM=0; addThis && clM<nbCats(); clM++)
		  {
		     for(int m=0; 
			 m<ds[clM].distinct() && ds[clM](m,0)<min0;
			 m++);
		     for(; addThis && m<ds[clM].distinct() && 
			    ds[clM](m,0)<=max0; 
			 m++)
			addThis=	(clM==clL && m==l) 
			   || (clM==clR && m==r) 
			   ||	!isInBox(ds[clL](l),ds[clR](r),
			   ds[clM](m));
		  }
	       }
	       if ( addThis )
	       {
		  list(ptr,0)=clL; list(ptr,1)=l; 
		  list(ptr,2)=clR; list(ptr,3)=r;
		  ptr++; 
	       }
	       else if ( trace>2 )
		  cout<< endl << "point " << ds[clL].label(l) 
		     << " in class " << clL
		     << " and point " << ds[clR].label(r) 
			<< " in class " << clR
			<< "suppressed by point-in-a-box";
	    }
	 }
   if ( trace>1 && ptr<list.m() )
      cout<< "\n*** Point-in-a-box method reduces the number of pairs from "
	  << list.m() << " to " << ptr << ".\n";
   return list.resize_(ptr,4);
}



tcT multiDS<T>& multiDS<T>::binarize(const cutPtsSet<T>& cps)
{
   int garbage;
   Matrix<char> oldMono(ds[0].monotone);
   for (int set=0; set<nbCats(); set++)
   {
      Matrix<T> bd(ds[set].distinct(),cps.nbCP());
      for (int row=0; row<bd.m(); row++)
	 for (int t=0; t<bd.n(); t++)
	    if(unknown(ds[set](row,cps.attribute(t))))
	       bd(row,t)=dont_know();
	    else bd(row,t) = cps.above(ds[set](row,cps.attribute(t)), t);
      ds[set] = dataSet<T>(bd,ds[set]);
      ds[set].checkMult(garbage).sort().reorganize();
      normalized=true;
   }
   if ( oldMono.n() )
   {
      Matrix<char> newMono(cps.nbCP());
      for (int cp=0; cp<cps.nbCP(); cp++)
	 newMono(cp) = oldMono(cps.attribute(cp));
      for (set=0; set<nbCats(); set++)
	 ds[set].monotone.resize(newMono);
   }
   return *this;
}



tcT multiDS<T>& multiDS<T>::binarize(const cutPtsSet<T>& cps, float confidence)
{
   int garbage;
   Matrix<char> oldMono(ds[0].monotone);
   for (int set=0; set<nbCats(); set++)
   {
      Matrix<T> bd(ds[set].distinct(),cps.nbCP());
      for (int row=0; row<bd.m(); row++)
	 for (int t=0; t<bd.n(); t++)
	 {
	    bd(row,t) = cps.compare(ds[set](row,cps.attribute(t)),t,confidence);
	    if (bd(row,t)==-1) bd(row,t)=dont_know();
	 }
      ds[set] = dataSet<T>(bd,ds[set]);
      ds[set].checkMult(garbage).sort().reorganize();
      normalized=true;
   }
   if ( oldMono.n() )
   {
      Matrix<char> newMono(cps.nbCP());
      for (int cp=0; cp<cps.nbCP(); cp++)
	 newMono(cp) = oldMono(cps.attribute(cp));
      for (set=0; set<nbCats(); set++)
	 ds[set].monotone.resize(newMono);
   }
   return *this;
}



const int percPerturbe=5;
float Perturbe(float x) { return x + x * percPerturbe * (2*randomR()-1)/100; }


tcT void multiDS<T>::mapCube (Matrix<float>& net, const int expectedDim/**=10**/)
{
   // define 2 to the power expectDim
   int twoToDim=1;
   for (int i=0; i<expectedDim; i++) twoToDim *= 2;

   // construct an index table for each point (without multiplicity)
   // as well as the mean vector in the input space
   Matrix<int> index(card(),2), ind(card(),2);
   Matrix<float> meanData(1,dim()); meanData=0.0;
   int pos=0;
   for (int set=0; set<nbCats(); set++)
      for (int pt=0; pt<ds[set].distinct(); pt++)
      {
	 for (int m=0; m<ds[set].multiplicity(pt); pos++, m++)
	 {
	    index(pos,0)=set;
	    index(pos,1)=pt;
	 }
	 for (int d=0; d<dim(); d++)
	    meanData(d) += float(ds[set](pt,d)) * ds[set].multiplicity(pt);
      }
   meanData /= card();
   
   // Define and initialize the network
   net.resize(twoToDim,dim());
   net = meanData;
   net = net.map(&Perturbe);

   // Learning rate parameter
   float eta0; cout<< "Give the initial learning rate 0:"; cin>> eta0;
   float decr0;cout<< "Give the decr factor for rate 0:"; cin>> decr0;
   float eta1; cout<< "Give the initial learning rate 1:"; cin>> eta1;
   float decr1;cout<< "Give the decr factor for rate 1:"; cin>> decr1;
   float eta2; cout<< "Give the initial learning rate 2:"; cin>> eta2;
   float decr2;cout<< "Give the decr factor for rate 2:"; cin>> decr2;
   float lastEta=0.1;
   
   // Main loop constructing the network
   Matrix<float> point(dim());
   int minDist,dist,winner;
   Matrix<float> clas(twoToDim); clas=0;
   int clasW=50; cout<<"Class weight:"; cin>> clasW;
#if _graphic
   const int xmin=0; const int xmax=40; const int ymin=0; const int ymax=40;
   const int offs=3; 
   const int xr1=xmin-offs/2; const int yr1=ymin-offs/2;
   const int xr2=xmax+offs/2; const int yr2=ymax+offs/2;
   const int xtext=xmin-offs+1; const int ytext=ymin-offs+1;
   const int sca=100;
   char symb[4]; symb[0]='o'; symb[1]='\0'; symb[2]='+'; symb[3]='\0';
   openpl(); erase(); 
   space(sca*(xmin-offs),sca*(ymin-offs),sca*(xmax+offs),sca*(ymax+offs)); 
   box(sca*xr1,sca*yr1,sca*xr2,sca*yr2); closepl();
#endif 
   do
   {
#if _graphic
      if ( dim()==2 )
      {
	 erase(); box(sca*xr1,sca*yr1,sca*xr2,sca*yr2);
	 for (int j=0; j<index.m(); j++)
	 {
	    move(sca*ds[index(j,0)](index(j,1),0),
		 sca*ds[index(j,0)](index(j,1),1));
	    ::label(&symb[2*index(j,0)]);
	 }
	 closepl();
	 for (i=0; i<twoToDim; i++)
	 {
	    circle(sca*net(i,0),sca*net(i,1),sca/10);
	    for (int bit=1; bit<twoToDim; bit<<=1)
	       if (!(bit & i))
		  line(sca*net(i,0),sca*net(i,1),
		       sca*net(i^bit,0),sca*net(i^bit,1));
	 }
	 closepl();
      }
#endif 
      ind=index;
      for (int nbpt=card(); nbpt>0; nbpt--)
      {
	 // pick one point at random
	 int temp = randomI(nbpt);
	 set = ind(temp,0); int pt = ind(temp,1);
	 ind.s(temp) = ind.s(nbpt-1);
	 for (int d=0; d<dim(); d++)
	    point(d) = float(ds[set](pt,d));
	 // determine the winner
	 minDist=10000;
	 for (i=0; i<twoToDim; i++)
	 {
	    dist=clasW*Sqr(2*set-1-clas(i));
	    for (d=0; d<dim(); d++) 
	      dist += Sqr(point(d)-net(i,d));
	    if (minDist > dist) { minDist=dist; winner=i; }
	 }
	 // modify the winner
	 net.s(winner) += eta0 * (point - net.s(winner));
	 clas(winner) += eta0 * (2*set-1 - clas(winner));
	 // modify the neighbours
	 int neighbour;
	 for (int bit=1; bit<twoToDim; bit<<=1)
	 {
	    neighbour = winner ^ bit;
	    net.s(neighbour) += eta1 * (point - net.s(neighbour));
	    clas(neighbour) += eta1 * (2*set-1 - clas(neighbour));
	    int neighbour2;
	    for (int bit2=bit<<1; bit2<twoToDim; bit2<<=1)
	    {
	       neighbour2 = neighbour ^ bit2;
	       net.s(neighbour2) += eta2 * (point - net.s(neighbour2));
	       clas(neighbour2) += eta2 * (2*set-1 - clas(neighbour2));
	    }
	 }

      }
      eta1*=decr1; eta2*=decr2;
   } while ((eta0*=decr0) > lastEta);
}



tcT multiDS<T>& multiDS<T>::binarize(const Matrix<float> net)
{
   int garbage;
   int oldDim=dim();
   for (int set=0; set<nbCats(); set++)
   {
      // Find the size of the binarized dataset
      for (int twoToDim=1, newDim=0; twoToDim!=net.m() && newDim<30; 
	   newDim++, twoToDim*=2);
      if ( newDim>=30 ) 
	 throw whereWhich(
	    "multiDS<T>& multiDS<T>::binarize(const Matrix<float> net)",
	    "Number of columns of net is not a power of 2");
      Matrix<T> bd(ds[set].distinct(),newDim);
      for (int row=0; row<bd.m(); row++)
      {
	 // Determine the winner
	 int minDist,dist,winner;
	 minDist=10000;
	 for (int i=0; i<twoToDim; i++)
	 {
	    dist=0;
	    for (int d=0; d<oldDim; d++)
	       dist+=(ds[set](row,d)-net(i,d))*(ds[set](row,d)-net(i,d));
	    if (minDist > dist) { minDist=dist; winner=i; }
	 }
	 // Binarize the point indexed by row
	 for (int d=0; d<newDim; d++, winner>>=1)
	    bd(row,d) = winner & 1;
      }
      ds[set] = dataSet<T>(bd,ds[set]);
      ds[set].checkMult(garbage).sort().reorganize();
      normalized=true;
   }
   return *this;
  
}


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


int patInd (int n, int d, binArray& wordP, const Matrix<int>& binomial)
{
   if ( n==d || !d ) return 0;
   int x1;
   for (x1=0; x1<wordP.n(); x1++) if ( wordP(x1) ) break;
   if ( d==1 ) return x1;
   return binomial(d,n) - binomial(d,n-x1) 
      + patInd(n-x1-1,d-1,wordP.s(0,1,x1+1),binomial);
}


void indPat (int n, int d, int index, binArray& wordP, 
	     const Matrix<int>& binomial)
{
   if ( !d ) wordP = false;
   else if ( d==1 ) { wordP = false; wordP.setOn(index); }
   else if ( n==d ) wordP = true;
   else
   {
      int nextN = n-1;
      while ( index >= binomial(d-1,nextN) ) index-=binomial(d-1,nextN--);
      if ( nextN < n-1 ) wordP.s(0,1,0,n-1-nextN)=false;
      wordP.setOn(n-1-nextN);
      indPat(nextN,d-1,index,wordP.s(0,1,n-nextN),binomial);
   }
}


void nextWord (binMatrix& w)
{
   int l=w.n()-1;
   if ( w(l) )
   {
      while ( l && w(l-1) ) l--;
      int nb1=w.n()-l;
      while ( l && !w(l-1) ) l--;
      w.s(0,1,l+1,nb1) = true;
      if ( l+1+nb1 < w.n() )
	 w.s(0,1,l+1+nb1) = false;
   }
   else while ( l && !w(l-1) ) l--;
   if ( !l ) throw whereWhich("void nextWord(binMatrix&)","No next word");
   w.setOn(l);
   w.setOff(l-1);
}


boolean coverOne(const Matrix<int>& rowInd, const binArray& table)
{
   unsigned char conjunction;
   for (int col=0; col<ceilDB(table.n()); col++)
   {
      conjunction = table.block(rowInd(0),col);
      for (int row=1; row<rowInd.n(); row++)
      {
	 conjunction &= table.block(rowInd(row),col);
	 if ( !conjunction ) break;
      }
      if ( conjunction ) return true;
   }
   return false;
}


boolean coverAtLeast(const Matrix<int>& rowInd, const binArray& table, 
		     const Array<int>& mult, const int minCover)
{
   char unsigned conjunction;
   int nbCover=0;
   for (int col=0; col<ceilDB(table.n()); col++)
   {
      conjunction = table.block(rowInd(0),col);
      for (int row=1; row<rowInd.n(); row++)
      {
	 conjunction &= table.block(rowInd(row),col);
	 if ( !conjunction ) break;
      }
      if ( conjunction )
      {
	 int stop=minimum(packetSize,table.n()-col*packetSize);
	 for (row=0; row<stop; row++)
	    if ( conjunction & p2[row] ) 
	       nbCover+=mult(packetSize*col+row);
      }
      if ( nbCover>=minCover ) return true;
   }
   return false;
}


int nbCover(const Matrix<int>& rowInd, const binArray& table, 
	    const Array<int>& mult)
{
   char unsigned conjunction;
   int res=0;
   for (int col=0; col<ceilDB(table.n()); col++)
   {
      conjunction = table.block(rowInd(0),col);
      for (int row=1; row<rowInd.n(); row++)
      {
	 conjunction &= table.block(rowInd(row),col);
	 if ( !conjunction ) break;
      }
      if ( conjunction )
      {
	 int stop=minimum(packetSize,table.n()-col*packetSize);
	 for (row=0; row<stop; row++)
	    if ( conjunction & p2[row] ) 
	       res+=mult(packetSize*col+row);
      }
   }
   return res;
}


void whoIsCovered(const Matrix<int>& rowInd, const binArray& table, 
		  Matrix<int>& coverage)
{
   binMatrix conjunction(table.n());
   coverage.resize(table.n());
   int nbC=0;
   conjunction = table.s(rowInd(0));
   for (int row=1; row<rowInd.n(); row++)
      conjunction &= table.s(rowInd(row));
   for (int i=0; i<conjunction.n(); i++)
      if ( conjunction(i) ) coverage(nbC++)=i;
   coverage.resize_(nbC);
}


boolean minDistAtLeast(const Matrix<int>& rowInd, const binArray& table, 
		       const int minDist, int toler)
{
   for (int col=0; col<table.n(); col++)
   {
      int dist=0;
      for (int row=0; row<rowInd.n(); row++)
	 if ( !table(rowInd(row),col) ) 
	    if ( ++dist >= minDist ) break;
      if ( dist<minDist )
	 if ( !toler )
	    return false;
	 else
	    toler--;
   }
   return true;
}


void update(
   const binMatrix& wordP, int wordI, int barsI, int n, int size, int maxSize,
   const binArray& sop, const binMatrix& pos, const binMatrix& neg,
   const Array<int>& sopMult, const Array<int>& posMult, 
   const Array<int> negMult, const int nbPos, const int nbNeg,
   Matrix<int>& posCover, int gap2Neg, float tolerance, binMatrix& might, 
   int& nbCandidates, int minPosCoverage, patterns& pList, int& nbPatterns)
{
   Matrix<int> indices(size);
   int ic=0; int bars1I=barsI;
   for (int pp=0; pp<n; pp++)
      if ( wordP(pp) )
      {
	 if ( bars1I % 2 ) 	indices(ic++)=pp+n;
	 else				indices(ic++)=pp;
	 bars1I/=2;
      }
   int nbCov;
   if ( (nbCov=nbCover(indices,sop,sopMult)) >= minPosCoverage )
   {
      int toler=(int)(float(nbCov)*tolerance*float(nbNeg)/float(nbPos));
      if ( coverAtLeast(indices,neg,negMult,toler+1) ||
	   (gap2Neg>1 && !minDistAtLeast(indices,neg,gap2Neg,toler)) )
      {
	 if ( size<maxSize ) might.setOn(wordI,barsI);
	 nbCandidates++;
      }
      else
      {
	 Matrix<int> covered;
	 whoIsCovered(indices,pos,covered);
	 for (int Iam=0; Iam<covered.n(); Iam++)
	    posCover(covered(Iam))++;
	    pList.add(n,indices,nbCover(indices,pos,posMult));
	 nbPatterns++;
      }
   }
}


tcT void multiDS<T>::patchPatterns
(patterns& patList, const Array<int>& posClass, const Array<int>& negClass, 
 const boolean orientation, int& puc) const
{
   int i,j,row,set,record;

   int posDistinct=0, max_dist=0; 
   int n = dim();
   for (i=0; i<posClass.n(); i++) 
   {
      posDistinct+=ds[posClass(i)].distinct();
      if(max_dist< ds[posClass(i)].distinct()) 
	 max_dist = ds[posClass(i)].distinct();
   }
   Matrix<int> pos_uncover(posClass.n(),max_dist);
   for (set=0; set<posClass.n(); set++)
      for (record=0; record<ds[posClass(set)].distinct(); record++)
      {
	 Array<T> thisData(ds[posClass(set)](record));
	 pos_uncover(set,record)=1; 
	 for (patList.first(); !patList.eol(); patList.next())
	    if ( patList.fire(thisData) ) 
	    { pos_uncover(set,record)=0; break; }
      }

   int cur_row=-1, cur_set=-1;
   for(set=0; set<posClass.n(); set++)
      for(row=0;row< ds[posClass(set)].distinct(); row++)
	 if (pos_uncover(set,row)==1) {cur_row=row; cur_set=set; break;}
   while (cur_row>=0)		
      // all this is for the positive patterns
   {
      Matrix<int> candidate(1,2*n,0);
      for(j=0;j<n;j++)				
	 // here we make the pattern to cover only the current uncovered point
	 if ( ds[posClass(cur_set)](cur_row,j)==1 )
	 {
	    if ( orientation ? monotone(j)!=negativeAttr
		             : monotone(j)!=positiveAttr )
	    candidate(j)=1;
	 }
	 else if ( ds[posClass(cur_set)](cur_row,j)==0 )
	 {
	    if ( orientation ? monotone(j)!=positiveAttr
		             : monotone(j)!=negativeAttr )
	    candidate(n+j)=1;
	 }
      
      int literal, maxMinDistance;

      // Check whether this pattern covers some negative points. In principle 
      // it should not, but if there are some monotonicity constraints, it can.
      boolean covers1Neg=false;
      for(set=0; set<negClass.n(); set++)
      {
	 for(int k=0;k<ds[negClass(set)].distinct();k++) 
	 {
	    int distance=0;
	    for(int c=0;c<n;c++)
	       if(ds[negClass(set)](k,c)==1 && candidate(n+c)==1) 
		  ++distance;
	       else if (ds[negClass(set)](k,c)==0 && candidate(c)==1) 
		  ++distance;
	    if ( !distance )
	    {
	       covers1Neg=true; break;
	    }
	 }
	 if ( covers1Neg ) break;
      }
      if ( covers1Neg )
      {
	 flog<<"\n### Data" << setw(5) << ds[posClass(cur_set)].label(cur_row)
	     << " cannot be covered by any pattern respecting monotonicity "
	     << "constraints (in void multiDS<T>::patchPattern) !\n";
	 // get this point out of the list of uncovered points
	 pos_uncover(cur_set,cur_row)=0;
      }
      else
      {
	 do
	 {
	    literal=0, maxMinDistance=0;
	    for(j=0;j<2*n; j++)		
	       // here we find the literal to be dropped next
	       if(candidate(j)) 
	       {
		  int minDistance=n+1;
		  candidate(j)=0;		
		  // here we momentarily drop the candidate literal
		  for(set=0; set<negClass.n(); set++)
		     for(int k=0;k<ds[negClass(set)].distinct();k++) 
		     {
			int distance=0;
			for(int c=0;c<n;c++)
			   if(ds[negClass(set)](k,c)==1 && candidate(n+c)==1) 
			      ++distance;
			   else if (ds[negClass(set)](k,c)==0 && candidate(c)==1) 
			      ++distance;
			if (distance < minDistance) {minDistance=distance;}
		     }
		  if (minDistance > maxMinDistance) 
	          {literal=j; maxMinDistance= minDistance;}
		  candidate(j)=1;		
		  // put the current literal back to try with the next
	       }
	    if (maxMinDistance>0) candidate(literal)=0;
	    // here we drop the chosen literal 
	 } while(maxMinDistance>0);		
	 // keep dropping literals while possible
	 
	 int size=0,next=0;	
	 // here we measure the size of the pattern to be included
	 for(i=0;i<n;i++) if(candidate(i)==1 || candidate(i+n)==1) ++size;
	 Matrix<int> index(size);
	 // here we put the pattern in the index matrix
	 for(i=0;i<2*n;i++) if(candidate(i)==1) index(next++)=i;
	 
	 int cover=0,fire;		
	 // here we find the coverage of the new pattern
	 for (set=0; set<posClass.n(); set++)
	    for (record=0; record<ds[posClass(set)].distinct(); record++)
	    {
	       
	       fire=1;
	       for(int c=0; c< dim(); c++)
		  if (ds[posClass(set)](record,c)!=1 && candidate(c)==1) 
	          {fire=0; break;}
		  else if (ds[posClass(set)](record,c)!=0 && candidate(n+c)==1) 
	          {fire=0; break;}
	       if(fire==1) 
	       { 	++cover;
	       // for(i=0;i<n;i++) cout << ds[posClass(set)](record,i) << " "; 
	       //cout << "\n";
	       }
	    }
	 
	 patList.add(n,index,cover);		
	 // here we add the new pattern to the list
	 puc=0;
	 for (set=0; set<posClass.n(); set++)
	    for (record=0; record<ds[posClass(set)].distinct(); record++)
	    {
	       Array<T> thisData(ds[posClass(set)](record)); ++puc;
	       for (patList.first(); !patList.eol(); patList.next())
		  if ( patList.fire(thisData) ) 
		  { pos_uncover(set,record)=0; --puc; break; }
	    }
      }
      cur_row=-1, cur_set=-1;
      for(set=0; set<posClass.n(); set++)
	 for(row=0;row< ds[posClass(set)].distinct(); row++)
	    if (pos_uncover(set,row)==1) {cur_row=row; cur_set=set; break;}
   }							
   // end of procedure for positive patterns

}
// End of function patchPatterns 


// Begin function printpatfile
tcT void multiDS<T>::printPatFile
(patterns& posList, patterns& negList, const Array<int>& posClass,
 const Array<int>& negClass, ostream& s, int getPat, int pos_neg, 
 Array<int>& stat) const
{
   int set, row, c_pos, e_pos, c_neg, e_neg, counter=0;
   int max_dist=0, i;
   int pos_count[2], neg_count[2];
   double total_pos, total_neg;
   pos_count[0]=pos_count[1]=neg_count[0]=neg_count[1]=0;
   for (i=0; i<posClass.n(); i++) 
   {
      if(max_dist< ds[posClass(i)].distinct()) 
	 max_dist = ds[posClass(i)].distinct();
   }
   Matrix<int> pos_error(posClass.n(),max_dist);
   max_dist=0;
   for (i=0; i<negClass.n(); i++) 
   {
      if(max_dist< ds[negClass(i)].distinct()) 
	 max_dist = ds[negClass(i)].distinct();
   }
   Matrix<int> neg_error(negClass.n(),max_dist);
   for(set=0;set<posClass.n(); set++) 
      // here we find out if a positive point is misclasified
      for (row=0; row< ds[posClass(set)].distinct(); row++)
      {
	 total_pos=total_neg=0;
	 Array<T> thisData(ds[posClass(set)](row));
	 pos_error(set,row)=0; 
	 pos_count[1]+=ds[posClass(set)].multiplicity(row);
	 for (posList.first(); !posList.eol(); posList.next())
	    if ( posList.fire(thisData) ) 
	    { total_pos += posList.weight(); }
	 for (negList.first(); !negList.eol(); negList.next())
	    if ( negList.fire(thisData) ) 
	    { total_neg += negList.weight(); }
	 if(total_pos <= total_neg) 
	 {
	    pos_error(set,row)=1;
	    pos_count[1]-=ds[posClass(set)].multiplicity(row);
	    pos_count[0]+=ds[posClass(set)].multiplicity(row);
	 }
      }

   for(set=0;set<negClass.n(); set++) 
      // here we find out if a negative point is misclasified
      for (row=0; row< ds[negClass(set)].distinct(); row++)
      {
	 total_pos=total_neg=0;
	 Array<T> thisData(ds[negClass(set)](row));
	 neg_error(set,row)=0; 
	 neg_count[1]+=ds[negClass(set)].multiplicity(row);
	 for (posList.first(); !posList.eol(); posList.next())
	    if ( posList.fire(thisData) ) 
	    { total_pos += posList.weight(); }
	 for (negList.first(); !negList.eol(); negList.next())
	    if ( negList.fire(thisData) ) 
	    { total_neg += negList.weight(); }
	 if(total_neg <= total_pos) 
	 {
	    neg_error(set,row)=1;
	    neg_count[1]-=ds[negClass(set)].multiplicity(row);
	    neg_count[0]+=ds[negClass(set)].multiplicity(row);
	 }
      }

   if(!getPat) 
   {
      s <<  "           Training                       Testing             Patterns  \n" 
	<<  "  c_pos  e_pos  c_neg  e_neg    c_pos  e_pos  c_neg  e_neg  ";
      if(pos_neg) 	s << "  positive\n\n"; else 	s << "  negative\n\n";
   }
   if (pos_neg) 
   {  
      stat(0,0)= pos_count[1]; stat(0,1)= pos_count[0];
      stat(0,2)= neg_count[1]; stat(0,3)= neg_count[0];
   }
   else 
   {
      stat(0,0)= neg_count[1]; stat(0,1)= neg_count[0];
      stat(0,2)= pos_count[1]; stat(0,3)= pos_count[0];
   }
   for (posList.first(); !posList.eol(); posList.next())
   {
      counter++;
      c_pos=e_pos=c_neg=e_neg=0;
      for(set=0; set< posClass.n(); set++)
	 for (row=0; row< ds[posClass(set)].distinct(); row++)
	 {
	    Array<T> thisData(ds[posClass(set)](row));
	    if(posList.fire(thisData) && pos_error(set,row)==0) 
	       c_pos+=ds[posClass(set)].multiplicity(row);
	    else if(posList.fire(thisData) && pos_error(set,row)==1) 
	       e_pos+=ds[posClass(set)].multiplicity(row);
	 }
      for(set=0; set< negClass.n(); set++)
	 for (row=0; row< ds[negClass(set)].distinct(); row++)
	 {
	    Array<T> thisData(ds[negClass(set)](row));
	    if(posList.fire(thisData) && neg_error(set,row)==0) 
	       c_neg+=ds[negClass(set)].multiplicity(row);
	    else if(posList.fire(thisData) && neg_error(set,row)==1) 
	       e_neg+=ds[negClass(set)].multiplicity(row);
	 }
      if(pos_neg)
      {
	 stat(counter,0)= c_pos; 
	 stat(counter,1)= e_pos;
	 stat(counter,2)= c_neg;
	 stat(counter,3)= e_neg;
      }
      else 
      {
	 stat(counter,0)= c_neg; 
	 stat(counter,1)= e_neg;
	 stat(counter,2)= c_pos;
	 stat(counter,3)= e_pos;
      }

   }



}
// End function printpatfile



tcT void multiDS<T>::smallPatterns
(  patterns& patList, int& nbNonCovered,
   const Array<int>& posClass, const Array<int>& negClass,
   const boolean orientation,
   int maxSize/**=4**/, int minTrue/**=1**/, int coverNeeded/**=10**/, 
   int gap2F/**=1**/, float tolerance/**=0.0**/  //, float rateLimit/**=0.3**/
) const
{
   chrono smallPatGen;
   int n = dim();
   int i,j,record;
   int posCard=0; int posDistinct=0; int negCard=0; int negDistinct=0;
   for (i=0; i<posClass.n(); i++) 
   {
      posCard+=ds[posClass(i)].card();
      posDistinct+=ds[posClass(i)].distinct();
   }
   for (i=0; i<negClass.n(); i++) 
   {
      negCard+=ds[negClass(i)].card();
      negDistinct+=ds[negClass(i)].distinct();
   }
   boolean reduceTrue = (3*coverNeeded < posCard);
   if ( minTrue<0 )
      minTrue = maximum(1,int(float(-minTrue)*float(posCard)/100.0));

   // construct two boolean matrices transposed and with the chosen order
	
   binMatrix pos(2*n,posDistinct,false);
   binMatrix sop(2*n,posDistinct,false);
   binMatrix neg(2*n,negDistinct,false);
   Matrix<int> posM(posDistinct);
   Matrix<int> sopM(posDistinct);
   Matrix<int> negM(negDistinct);
   for (i=0; i<n; i++)
   {
      int col=0;
      for (int set=0; set<negClass.n(); set++)
	 for (record=0; record<ds[negClass(set)].distinct(); record++, col++)
	 {
	    if ( ds[negClass(set)](record,i)!=0 ) neg.setOn(i,col);
	    if ( ds[negClass(set)](record,i)!=1 ) neg.setOn(n+i,col);
	    negM(col)=ds[negClass(set)].multiplicity(record);
	 }
      col=0;
      for (set=0; set<posClass.n(); set++)
	 for (record=0; record<ds[posClass(set)].distinct(); record++, col++)
	 {
	    if ( ds[posClass(set)](record,i)==1 ) pos.setOn(i,col);
	    if ( ds[posClass(set)](record,i)==0 ) pos.setOn(n+i,col);
	    posM(col)=ds[posClass(set)].multiplicity(record);
	 }
   }
   sop = pos; sopM = posM;
   int nbSop = sop.n();
   Matrix<int> posCover(1,posDistinct,0);
   Matrix<int> sop2pos(posDistinct);
   for (i=0; i<posDistinct; i++) sop2pos(i)=i;


   // Some useful tables
   Matrix<int> binomial(maxSize+1,n+1,1);
   Matrix<int> power2(maxSize+1); power2(0)=1;
   for (i=1; i<=maxSize; i++)
   {
      power2(i) = 2*power2(i-1);
      for (j=i+1; j<=n; j++)
	 binomial(i,j)=binomial(i-1,j-1)+binomial(i,j-1);
   }


   // Declaration of the global variables containing the patterns
   // And the potential patterns

   binMatrix mightbcp[2];
   int current,previous;
   if ( maxSize < 2 )
   {
      mightbcp[0].resize(1,1);
      mightbcp[1].resize(1,1);
      current = 0;
   }
   else if ( 2*maxSize -4 < n )
   {
      mightbcp[0].resize(binomial(maxSize-2,n),power2(maxSize-2));
      mightbcp[1].resize(binomial(maxSize-1,n),power2(maxSize-1));
      current = maxSize%2;
   }
   else 
   {
      mightbcp[0].resize(binomial(n/2,n),power2(maxSize-2));
      mightbcp[1].resize(binomial(n/2,n),power2(maxSize-1));
      current = maxSize%2;
   }
   mightbcp[current].s(0,1,0,1) = true;

   patList.empty();
   float candidateRate = 1.0;

   if ( trace>2 )
   {
      cout << "\n  time #Po dg #patt #candid    #terms |";
      for (i=0; i<21; i++) cout << setw(3) << i; cout << "  >\n";
      cout.precision(1); cout.setf(ios::showpoint | ios::fixed);
      cout<< setw(6) << smallPatGen.time() << setw(4) << nbSop << setw(3) 
	 << 0 << setw(6) << 0 << setw(8) << 1 << setw(10) << 1 << " |\n";
   }

   // Main loop on the size of the terms

   for (int size=1; size<=maxSize; size++)
   {
      previous = current;
      current  = 1-current;
      if ( size < maxSize )
      {
	 // initialization of the table with patterns of size SIZE
	 // This table is initialized empty
	 mightbcp[current].s(0,binomial(size,n),0,power2(size)) = false;
      }
		
      // declaration of some working variables
      int nbPatterns=0; int nbCandidates=0;

      if ( size > 1 && false ) // candidateRate <= rateLimit )
      {
	 boolean foundOne = false;
	 binMatrix wordP(1,n);
	 binMatrix word1P(1,n);
	 int fromVar;

	 int nbBlocks = ( size < 5 ) ? 1 : power2(size-4);

	 // loop on the word (sequence of variables) of size SIZE-1
	 for (int wordI=0; wordI<binomial(size-1,n); wordI++)
	 {
	    foundOne = false;

	    // loop on the possible bars, by blocks of 8
	    for (int barsB=0; barsB<nbBlocks; barsB++)
	       if ( mightbcp[previous].block(wordI,barsB) )
	       {
		  if ( !foundOne ) // first found for this word
		  {
		     foundOne=true;
		     indPat(n,size-1,wordI,wordP,binomial); // compute wordP
		     if ( wordI != patInd(n,size-1,wordP,binomial) )
			throw whereWhich(
			   "void multiDS<T>::smallPatterns(...)",
			   "patInd(indPat(x))!=x");
		     if ( wordP(n-1) ) goto nextWordI;
		     if ( size==1 ) fromVar=0;
		     else
		     {
			fromVar=n-1;
			while ( !wordP(fromVar-1) ) fromVar--;
		     }
		  }
		  // at this point, foundOne=true
		  // and wordP and fromVar are defined
		  // loop on the possible bars within block barsB
		  int stopAt=minimum(packetSize*(barsB+1),power2(size-1));
		  for (int barsI=packetSize*barsB; barsI<stopAt; barsI++)
		     if ( mightbcp[previous](wordI,barsI) )
		     {
			// loop on all possible variables that can be added
			for (int newVar=fromVar; newVar<n; newVar++)
			{
			   // Check whether all its subterms are in 
			   // mightbcp[previous]
			   word1P = wordP;	word1P.setOn(newVar);
			   int rank=1;
			   int bars1I, rem;
			   boolean isCandidPos= 
			      (orientation ? monotone(newVar)!=negativeAttr
					   : monotone(newVar)!=positiveAttr);
			   boolean isCandidNeg=
			      (orientation ? monotone(newVar)!=positiveAttr
					   : monotone(newVar)!=negativeAttr);
			   int subInd;
			   const int barOS = power2(size-2);
			   // loop on all the subterms of the current word
			   for (int sub=0; sub<fromVar && 
				   (isCandidPos || isCandidNeg); sub++)
			      if ( word1P(sub) )
			      {
				 word1P.setOff(sub);
				 rem  = barsI % rank;
				 bars1I = rem + ((barsI-barsI%(rank*=2))/2);
				 subInd = patInd(n,size-1,word1P,binomial);
				 isCandidPos = isCandidPos && 
				    mightbcp[previous](subInd,bars1I);
				 isCandidNeg = isCandidNeg && 
				    mightbcp[previous](subInd,bars1I+barOS);
				 word1P.setOn(sub);
			      }

			   if ( isCandidPos || isCandidNeg )
			   {
			      int word1I = patInd(n,size,word1P,binomial);
			      if ( isCandidPos )
				 update(word1P,word1I,barsI,
					n,size,maxSize,
					sop.s(0,2*n,0,nbSop),pos,neg,sopM,posM,
					negM,posCard,negCard,posCover,gap2F,
					tolerance,mightbcp[current],
					nbCandidates,minTrue,
					patList,nbPatterns);
			      if ( isCandidNeg )
				 update(word1P,word1I,barsI+power2(size-1),
					n,size,maxSize,
					sop.s(0,2*n,0,nbSop),pos,neg,sopM,posM,
					negM,posCard,negCard,posCover,gap2F,
					tolerance,mightbcp[current],
					nbCandidates,minTrue,
					patList,nbPatterns);
			   }

			} // for newVar

		     } // if in for barsI
	       } // if in for basrB
	    nextWordI : ;
	 } // for wordI
				
      } // if ( candidateRate <= rateLimit )
      else // ( candidateRate > rateLimit )
      {
	 binMatrix wordP(1,n,false);
	 wordP.s(0,1,0,size) = true;
	 binMatrix word1P(1,n,false);
	 Matrix<int> subWordIndex(1,n);

	 // loop on the word (sequence of variables) of size SIZE
	 for (int wordI=0; wordI<binomial(size,n); wordI++)
	 {
	    if ( wordI ) nextWord(wordP);

	    subWordIndex = 0;

	    // loop on all possible combinations of negations
	    for (int barsI=0; barsI<power2(size); barsI++)
	    {
	       // Check whether the monotonicity constraints are fulfilled
	       if ( ds[0].monotone.n() ) // there are some monot. constraints
	       {
		  boolean monotonicityOK=true;
		  int barsA = barsI; int negated;
		  for (int v=0; v<wordP.n(); v++)
		     if ( wordP(v) )
		     {
			negated=barsA%2; barsA>>=1;
			if ( orientation ?
			       ( (  negated && monotone(v)==positiveAttr ) ||
				 ( !negated && monotone(v)==negativeAttr )
			       ) :
			       ( ( !negated && monotone(v)==positiveAttr ) ||
				 (  negated && monotone(v)==negativeAttr )
			       )
			   )
			{ monotonicityOK=false; break; }
		     }
		  if ( !monotonicityOK ) continue;
	       }
		  
	       // Check whether all its subterms are in mightbcp[previous]
	       word1P = wordP;
	       int rank=1;
	       int bars1I, rem;
	       boolean isCandid=true;
	       // loop on all the subterms of the current word
	       for (int sub=0; sub<n; sub++)
		  if ( word1P(sub) )
		  {
		     word1P.setOff(sub);
		     rem  = barsI % rank;
		     bars1I = rem + ((barsI-(barsI%(rank*=2)))/2);
		     if ( !subWordIndex(sub) )
			subWordIndex(sub) = 
			patInd(n,size-1,word1P,binomial);
		     if ( !mightbcp[previous](subWordIndex(sub),bars1I) )
		     {
			isCandid=false;
			break;
		     }
		     word1P.setOn(sub);
		  }
	       
	       if ( isCandid )
		  // Check whether some true points or false points 
		  // are covered
		  update(wordP,wordI,barsI,n,size,maxSize,
			 sop.s(0,2*n,0,nbSop),pos,neg,sopM,posM,negM,
			 posCard,negCard,
			 posCover,gap2F,tolerance,mightbcp[current],
			 nbCandidates,minTrue,patList,nbPatterns);
	       
	    } // for barsI
	 } // for wordI
      } // if ( candidateRate >  0.5 )

      if ( trace>2 )
      {
	 cout<< setw(6) << smallPatGen.time() << setw(4) << nbSop << setw(3) 
	    << size << setw(6) << nbPatterns << setw(8) << nbCandidates 
	    << setw(10) << power2(size)*binomial(size,n) << " |";
	 Matrix<int> distr(1,22,0);
	 for (i=0; i<posCover.n(); i++)
	    if (posCover(i)<21)
	       distr(posCover(i))++;
	    else
	       distr(21)++;
	       for (i=0; i<22; i++) 
		  if (distr(i)) 
		     cout << setw(3) << distr(i); 
		  else
		     cout << "   ";
	 cout << endl;
      }

      if ( !nbCandidates ) goto terminate;
		
      if ( reduceTrue && size < maxSize )
      {
	 // Check whether some points are enough covered
	 int prevNbPos = nbSop;
	 for (i=0; i<nbSop; i++)
	    if ( posCover(sop2pos(i)) >= coverNeeded )
	    {
	       sop.s(0,sop.m(),i,1) = sop.s(0,sop.m(),nbSop-1,1);
	       sopM(i) = sopM(nbSop-1);
	       sop2pos(i--) = sop2pos(--nbSop);
	    }
	 if ( !nbSop ) goto terminate;
	 if ( prevNbPos > nbSop )
	 {
	    if ( trace>2 )
	       cout<< setw(6) << smallPatGen.time() 
		  << setw(4) << nbSop << flush;
				
	    for (int wordI=0; wordI<binomial(size,n); wordI++)
	       for (int barsI=0; barsI<power2(size); barsI++)
		  if ( mightbcp[current](wordI,barsI) )
		  {
		     binMatrix wordP(n);
		     indPat(n,size,wordI,wordP,binomial);
		     Matrix<int> indices(size);
		     int ic=0; int bars1I=barsI;
		     for (int pp=0; pp<n; pp++)
			if ( wordP(pp) )
			{
			   if ( bars1I % 2 ) 	indices(ic++)=pp+n;
			   else				indices(ic++)=pp;
			   bars1I/=2;
			}
		     if ( !coverAtLeast(indices,sop.s(0,2*n,0,nbSop),
					sopM,minTrue) )
		     {
			mightbcp[current].setOff(wordI,barsI);
			nbCandidates--;
		     }
		  }

	    if ( trace>2 )
	       cout<< setw(17) << nbCandidates << endl;

	 } // if ( prevNbPos > nbSop )

      } // if ( reduceTrue && size<maxSize )

      candidateRate=float(nbCandidates)/float(power2(size)*binomial(size,n));
	
   } // for size

   terminate :
      nbNonCovered=0;
   for (i=0; i<posCover.n(); i++) if (!posCover(i)) nbNonCovered+=posM(i);
   patList.do_sort();
   //mightbcp[0].resize();
   //mightbcp[1].resize();
}


tcT void multiDS<T>::weightPatterns
  (patterns& patList, 
   const Array<int>& posClass, 
   const Array<int>& negClass,
   int weighting/**=1**/, 
   boolean normalize/**=true**/, 
   boolean delZWP/**=true**/) const
{
   switch(weighting)
   {
    case 0:
    for (patList.first(); !patList.eol(); patList.next())
       patList.weight()=1.0;
    break;

    case 1:
    for (patList.first(); !patList.eol(); patList.next())
       patList.weight()=float(patList.coverage());
    break;

    case 2:
    case 3:
    {
       int power2[10]; power2[0]=1;
       for (int i=1; i<10; i++) power2[i]=2*power2[i-1];
       if ( weighting==2 )
	  for (patList.first(); !patList.eol(); patList.next())
	     patList.weight()=
		float(power2[patList.size()]*patList.coverage());
       else // ( weighting==3 )
	  for (patList.first(); !patList.eol(); patList.next())
	     patList.weight()=float(power2[10-patList.size()]);
    }
    break;

    case 6:
    for (patList.first(); !patList.eol(); patList.next())
       patList.weight()=Sqr(float(patList.coverage()));
    break;

    case 7:
    // here the weight is the cube of the coverage Mario
    for (patList.first(); !patList.eol(); patList.next())
       patList.weight()=
	  Sqr(float(patList.coverage()))*(float(patList.coverage()));
    break;

    case 8:
    // here the weight is the exponentiation of the coverage. Mario
    for (patList.first(); !patList.eol(); patList.next())
    { 	
       float base;	// base for expontiation.  
       int i;
       base=1.0;	// value is arbritrary, small enough to prevent 
       // overflow (hopefully)
       for(i=1;i<=patList.coverage();i++) base*=1.2;  
       // so the weight is 1.2^coverage 
       patList.weight()=base;
    }
    break;

    case 4:
    {
       // count the number of points
       int nbPoints=0;
       for (int set=0; set<posClass.n(); set++)
       {
	  nbPoints+=ds[posClass(set)].distinct();
       }
       binMatrix covered(1,nbPoints,false);
       nbPoints=0; int finger=0;
       for (set=0; set<posClass.n(); set++)
	  for (int pt=0; pt<ds[posClass(set)].distinct(); pt++,finger++)
	  {
	     Array<T> thisData(ds[posClass(set)](pt));
	     for (patList.first(); !patList.eol(); patList.next())
		if ( patList.fire(thisData) )
		{
		   covered.setOn(finger);
		   nbPoints++;
		   break;
		}
	  }

       // construct the LP
       Matrix<int> c(1, patList.nbPatterns()+1, 0);
       c(patList.nbPatterns())=+1;
       Matrix<int> b(1, nbPoints+1, 0);
       b(nbPoints)=+1;
       Matrix<char> R(1, nbPoints+1, 'G');
       R(nbPoints)='E';
       Matrix<int> A(nbPoints+1, patList.nbPatterns()+1, 0);
       A.s(nbPoints,1,0,patList.nbPatterns())=+1;
       A.s(0,nbPoints,patList.nbPatterns(),1)=-1;
       finger=0; int row=0;
       for (set=0; set<posClass.n(); set++)
	  for (int pt=0; pt<ds[posClass(set)].distinct(); pt++, finger++)
	     if ( covered(finger) )
	     {
		Array<T> thisData(ds[posClass(set)](pt));
		int col=0;
		for (patList.first(); !patList.eol(); patList.next(), col++)
		   if ( patList.fire(thisData) )
		      A(row,col)=1;
		row++;
	     }

       // create and solve the LP
       linearProgram myLP (MAXIMIZE, c, A, R, b);
      
       double optValue;
       Matrix<double> optPoint(patList.nbPatterns()+1);
       int status = myLP.solve(optValue,optPoint);
       if ( status ) 
       {
	  flog<< "\n### lp status is " << status << "!!!\n" << flush; 
	  return;
       }
       int col=0;
       for (patList.first(); !patList.eol(); patList.next())
	  patList.weight()=optPoint(col++);

       normalize=false; // it is already normalized
    }
    break;
   } // switch (weighting)

   if ( delZWP )
   {
      float totWeight=0.0;
      for (patList.first(); !patList.eol(); patList.next())
	 totWeight+=patList.weight();
      float thr = (totWeight/float(patList.nbPatterns()))/10000.0;
      patList.first();
      while ( !patList.eol() )
	 if ( patList.weight() < thr )
	    patList.delPattern();
	 else
	    patList.next();
   }

   if ( normalize ) // normalize weights so that sum is 1.0
   {
      float totWeight=0.0;
      for (patList.first(); !patList.eol(); patList.next())
	 totWeight+=patList.weight();
      for (patList.first(); !patList.eol(); patList.next())
	 patList.weight()/=totWeight;
   }

   if ( trace>2 )
   {
      float mean=0.0; float min=1000.0; float max=-1000.0;
      for (patList.first(); !patList.eol(); patList.next())
      {
	 mean+=patList.weight(); 
	 minimize(min,patList.weight());
	 maximize(max,patList.weight());
      }
      mean /= float(patList.nbPatterns());
      cout<< "mean=" << setw(6) << setprecision(2) << mean 
	  << ",   span=" << setw(6) << (mean==0.0 ? 0.0 : (max-min)/mean)
	  << "   for weights of patterns covering class " << posClass;
   }
}


tcT void multiDS<T>::balanceWeightPatterns
  (patterns& posList, 
   patterns& negList, 
   const Array<int>& posClass, 
   const Array<int>& negClass,
   int method/**=1**/) const
{  // First, associate to each point (positive and negative) their multipli-
   // city as well as a triplet a = (f-, f+, s), where f- (resp. f+) is the
   // result of the pseudo-boolean function f+ (resp. f-), and s is +1 for
   // the positive points and -1 for the negative points.
   int distinctPos=0; int distinctNeg=0; int cardPos=0;
   for (int set=0; set<posClass.n(); set++)
   {
      distinctPos+=ds[posClass(set)].distinct();
      cardPos+=ds[posClass(set)].card();
   }
   for (set=0; set<negClass.n(); set++)
      distinctNeg+=ds[negClass(set)].distinct();
   Matrix<float> a(distinctPos+distinctNeg,3); a = 0.0;
   Matrix<int> m(1,distinctPos+distinctNeg);
   int point=0;
   for (set=0; set<negClass.n(); set++)
      for (int example=0; example<ds[negClass(set)].distinct(); 
	   example++, point++)
      {
	 m(point) = ds[negClass(set)].multiplicity(example);
	 Array<T> thisData(ds[negClass(set)](example));
	 for (negList.first(); !negList.eol(); negList.next())
	    if ( negList.fire(thisData) )
	       a(point,0)-=negList.weight();
	 for (posList.first(); !posList.eol(); posList.next())
	    if ( posList.fire(thisData) )
	       a(point,1)+=posList.weight();
	 a(point,2)=-1;
      }
   for (set=0; set<posClass.n(); set++)
      for (int example=0; example<ds[posClass(set)].distinct(); 
	   example++, point++)
      {
	 m(point) = ds[posClass(set)].multiplicity(example);
	 Array<T> thisData(ds[posClass(set)](example));
	 for (negList.first(); !negList.eol(); negList.next())
	    if ( negList.fire(thisData) )
	       a(point,0)-=negList.weight();
	 for (posList.first(); !posList.eol(); posList.next())
	    if ( posList.fire(thisData) )
	       a(point,1)+=posList.weight();
	 a(point,2)=1;
      }

   // Try to find a triplet x = (x+, x-, x0) such that a.x > 0 for all points
   // a. This is done by minimizing the criterion of the perceptron, using a
   // full gradient descent.
   Matrix<float> x(1,3); x(0)=-1.0; x(1)=1.0; x(2)=0.0; x /= sqrt((x % x)(0)); 
   Matrix<float> bestx(x); 
   int nbErrors=1;
   float objVal; float optVal=1E10;
   Matrix<float> grad(1,3);
   float step;   
   for (step=0.2; step>0.00001; step*=0.99)
   {
      nbErrors=0; grad=0.0;
      for (int pt=0; pt<a.m(); pt++)
	 if ( x % a.s(pt) <= 0.0 ) 
	 { nbErrors+=m(pt); grad += (float(m(pt)) * a.s(pt)); }
      objVal = -((x % grad)(0));
      if ( objVal < optVal ) { optVal = objVal; bestx = x; }
      if ( !nbErrors ) break;
      //grad /= sqrt(grad % grad);
      x += (step * grad); x /= sqrt((x % x)(0));
   }
   x = bestx;

	// if linearly separable, try to maximize the stability, still using a
	// gradient descent algorithm
   if ( !nbErrors )
   {
      optVal=0.0;
      float min;
      for (step=0.1; step>0.00001; step*=0.99)
      {
	 min=1.1*((x % a).min()(0)); grad=0.0;
	 for (int pt=0; pt<a.m(); pt++)
	    if ( x % a.s(pt) <= min ) grad += (float(m(pt)) * a.s(pt));
	 if ( min > optVal ) { optVal = min; bestx = x; }
	 //grad /= sqrt(grad % grad);
	 x += (step * grad); x /= sqrt((x % x)(0));
      }
   }
   x = bestx;

	// translate the solution x into an appropriate modification of the
	// weights of the patterns
   if ( x(0) >= 0.0)
      throw whereWhich(
	 "void multiDS<T>::balanceWeightPatterns(...) const",
	 "Strang solution x(0)>0");
   x /= -x(0);
   if ( x(1) != 1.0 )
      for (posList.first(); !posList.eol(); posList.next())
	 posList.weight() *= x(1);
   if ( x(2) != 0.0 )
   {
      Matrix<int> emptyPattern(0);
      posList.add(dim(),emptyPattern,cardPos,x(2));
   }							  
}


tcT	void multiDS<T>::cleanPatterns (patterns& patList, 
				      const Array<int>& posClass) const
{
   int nbDistinct=0;
   for (int set=0; set<posClass.n(); set++) 
      nbDistinct+=ds[posClass(set)].distinct();
   binMatrix A(patList.nbPatterns(),nbDistinct,false);
   int col=0; // point
   for (set=0; set<posClass.n(); set++)
      for (int point=0; point<ds[posClass(set)].distinct(); point++)
      {
	 boolean pointCovered=false;
	 Array<T> thisData(ds[posClass(set)](point));
	 int row=0;
	 for (patList.first(); !patList.eol(); patList.next(), row++)
	    if ( patList.fire(thisData) ) 
	    { A.setOn(row,col); pointCovered=true; }
	 if ( pointCovered ) col++;
      }
   A.resize_(A.m(),col);

   Matrix<int> info(patList.nbPatterns(),2); int i;
   for (i=0, patList.first(); i<patList.nbPatterns(); i++, patList.next()) 
   { info(i,0)=i; info(i,1)=patList.size(); } 
   int nbDel=0;
   Matrix<int> deleted(patList.nbPatterns());

   int nbPat=A.m(); int comp;
   for (int upper=0; upper<nbPat-1; upper++)
      for (int lower=upper+1; lower<nbPat; lower++)
	 if ( (comp=(A.s(upper) %= A.s(lower)))!=2 )
	    // the 2 rows are comparable
	    if ( ( comp+=tcomp(info(lower,1),info(upper,1)) ) > 0 )
	       // comp=a+b, a= 1 if A.s(upper) included in A.s(lower)
	       //			 a= 0 if A.s(upper) == A.s(lower)
	       //			 a=-1 if A.s(upper) contains A.s(lower)
	       //			 b= 1 if degree upper >  degree lower
	       //			 b= 0 if degree upper == degree lower
	       //			 b=-1 if degree upper <  degree lower
	       // Thus, if comp==+1,+2, upper is deleted, if comp==-1,-2,
	       // lower is deleted, while if comp==0, nobody is deleted.
	    {
	       A.s(upper) = A.s(--nbPat);
	       deleted(nbDel++)=info(upper,0);
	       info.s(upper--) = info.s(nbPat);
	       lower = nbPat;
	    }
	    else if ( comp < 0 )
	    {
	       A.s(lower) = A.s(--nbPat);
	       deleted(nbDel++)=info(lower,0);
	       info.s(lower--) = info.s(nbPat);
	    }

   if ( nbDel )
   {
      deleted.resize_(nbDel);
      deleted.sort();
      // delete the furthest in the list first, since delPattern
      // place the last element of the list at the place of the one
      // deleted.
      for (int ptr=nbDel-1; ptr>=0; ptr--) 
	 patList.delPattern(deleted(ptr));

      patList.do_sort();
   }
}


tcT	void multiDS<T>::reducePatterns ( patterns& patList, 
					const Array<int>& posClass, 
					int minCover/**=1**/) const
{
   int nbDistinct=0;
   for (int set=0; set<posClass.n(); set++) 
      nbDistinct+=ds[posClass(set)].distinct();
   binMatrix A(nbDistinct,patList.nbPatterns(),false);
   int minCoverage=patList.nbPatterns();
   int row=0; // one for each point
   for (set=0; set<posClass.n(); set++)
      for (int point=0; point<ds[posClass(set)].distinct(); point++)
      {
	 Array<T> thisData(ds[posClass(set)](point));
	 int col=0; int nbCover=0;
	 for (patList.first(); !patList.eol(); patList.next(), col++)
	    if ( patList.fire(thisData) ) 
	    { A.setOn(row,col); nbCover++; }
	 if ( nbCover ) 
	 {
	    row++;
	    minimize(minCoverage,nbCover);
	 }
	 if ( nbCover < minCover )
	    flog<< "\n### data " << setw(3) << ds[posClass(set)].label(point) 
		<< " is covered by " << nbCover << " pattern" 
		<< (nbCover ? (nbCover==1 ? " only" : "s only") : "") 
		<< endl << flush;
      }
   A.resize_(row,A.n());

   setCover sc(A);
   sc.greedy(flog,1,-minCover,false);
   if ( minCover==1 ) sc.reduce();	

   Matrix<int> sol(sc.solution());

   if ( sol.n() < A.n() )
   {
      sol.sort();
      // delete the furthest missing from the list first, since delPattern
      // place the last element of the list at the place of the one
      // deleted.
      int pos=sol.n()-1;
      for (int ptr=A.n()-1; ptr>=0; ptr--)
	 if (pos>=0 && ptr==sol(pos))
	    pos--;
	 else
	    patList.delPattern(ptr);
      patList.do_sort();
   }
}


tcT int multiDS<T>::nbErrors( patterns& posList, patterns& negList,
			      float errorTolerance, int& wronglyOn,
			      int& wronglyOff, int& undecidable,
			      char* fileName/**=NULL**/) const
{
   int res=0; wronglyOn=0; wronglyOff=0; undecidable=0;
   ofstream out, out_err;
   if ( fileName )
   {
      out.open(fileName);
      out.setf(ios::fixed);
      out.precision(5);
      char fne[80]; sprintf(fne,"%s.error",fileName);
      out_err.open(fne);
      out_err.setf(ios::fixed);
      out_err.precision(5);
   }
   for (int set=0; set<2; set++)
      for (int example=0; example<ds[set].distinct(); example++)
      {
	 float totalPos=0.0; float totalNeg=0.0;
	 Array<T> thisData(ds[set](example));
	 // evaluate the Pseudo-boolean function
	 for (posList.first(); !posList.eol(); posList.next())
	    if ( posList.fire(thisData) ) totalPos+=posList.weight();
	 for (negList.first(); !negList.eol(); negList.next())
	    if ( negList.fire(thisData) ) totalNeg+=negList.weight();
	 if ( fileName )
	    out	<< setw(1) << set << setw(4) << ds[set].label(example)
	       << setw(4) << ds[set].multiplicity(example) 
		  << setw(10) << totalPos
		  << setw(10) << totalNeg 
		  << setw(10) << totalPos-totalNeg
		  << endl;
	 if ( (  set && totalPos-totalNeg<-errorTolerance ) || 
	      ( !set && totalPos-totalNeg>+errorTolerance ) ) 
	 {
	    res+=ds[set].multiplicity(example);
	    if (set)
	       wronglyOff+=ds[set].multiplicity(example);
	    else
	       wronglyOn +=ds[set].multiplicity(example);
	    if ( fileName )
	    {
	       out_err	<< "\n\npoint " << ds[set].label(example)
		  << ", from class " << set
		  << setw(10) << totalPos
		  << setw(10) << totalNeg 
		  << setw(10) << totalPos-totalNeg;
	       int nbPinL=0;
	       for (posList.first(); !posList.eol(); posList.next())
		  if ( posList.fire(thisData) ) 
		  {
		     if (!nbPinL) 
			out_err << "\npositive firing patterns:\n";
		     else
			out_err << ",  ";
		     if (!(++nbPinL % 5)) out_err << endl;
		     out_err << "(";
		     posList.printPattern(out_err);
		     out_err << ")";
		  }
	       nbPinL=0;
	       for (negList.first(); !negList.eol(); negList.next())
		  if ( negList.fire(thisData) ) 
		  {
		     if (!nbPinL) 
			out_err << "\nnegative firing patterns:\n";
		     else
			out_err << ",  ";
		     if (!(++nbPinL % 5)) out_err << endl;
		     out_err << "(";
		     negList.printPattern(out_err);
		     out_err << ")";
		  }
	    }					
	 }
	 else if ( absolute(totalPos-totalNeg)<=errorTolerance )
	 {
	    undecidable+=ds[set].multiplicity(example);
	    res+=ds[set].multiplicity(example);
	 }
      }
   return res;
}



tcT int multiDS<T>::boundOnErrors(const multiDS<T>& basis) const
// In this form, this routine works only if basis is consistant.
// In case of inconsistant basis its a mess to decide what will be
// a mistake.
// We also assume that multiplicities are checked in basis and in *this.
{
   int res=0;
   if ( !basis.isConsistant() ) return -1;
   binMatrix error(1,distinct(),false);
   int cumul=0;
   for (int set=0; set<nbCats(); set++)
   {
      for (int ind=0; ind<ds[set].distinct(); ind++)
	 for (int setB=0; setB<basis.nbCats(); setB++)
	    if ( setB!=set )
	       for (int indB=0; indB<basis.ds[setB].distinct(); indB++)
		  if ( identical(ds[set](ind),basis.ds[setB](indB)) )
		  {
		     error.setOn(cumul+ind);
		     res+=ds[set].multiplicity(ind); 
		     setB=basis.nbCats(); break;
		  }
      cumul+=ds[set].distinct();
   }
   cumul=0;
   for (set=0; set<nbCats()-1; set++)
   {
      for (int ind=0; ind<ds[set].distinct(); ind++)
	 if ( !error(cumul+ind) )
	 {
	    int maxi = ds[set].multiplicity(ind);
	    res+=maxi;
	    int cumul1=cumul+ds[set].distinct();
	    for (int set1=set+1; set1<nbCats(); set1++)
	    {
	       for (int ind1=0; ind1<ds[set1].distinct(); ind1++)
		  if ( !error(cumul1+ind1) && 
		       ds[set](ind)==ds[set1](ind1) )
		  {
		     maxi = maximum(ds[set1].multiplicity(ind1),maxi);
		     res+=ds[set1].multiplicity(ind1);
		     break;
		  }
	       cumul1+=ds[set1].distinct();
	    }
	    res-=maxi;							
	 }
      cumul+=ds[set].distinct();
   }
   return res;
}


tcT dataSet<T> multiDS<T>::merge(const boolean class_first/**=true**/)
{
   Matrix <T> 	d(distinct(),dim()+1);// data set matrix
   Matrix <int> m(distinct());	      // multiplicity vector
   Matrix <int>	l(distinct());	      // label vector

   int position=0;
   for(int set=0;set< nbCats();position+=ds[set++].distinct())
   {
      d.s(position,ds[set].distinct(),class_first,dim())= 
	 ds[set].matrix();
      if(class_first) d.s(position,ds[set].distinct(),0,1)= (T) (set);
      else d.s(position,ds[set].distinct(),dim(),1)= (T) (set);
      for(int i=0; i<ds[set].distinct();i++)
      {
	 m(position+i) = ds[set].multiplicity(i);
	 l(position+i) = ds[set].label(i);
      }
   }
   return(dataSet<T>(d,m,l,ds[0].monotone));
}


tcT multiDS<T>& multiDS<T>::separate (const dataSet<T>& from)
{
#ifdef _CHECK_INDEX
   if ( from.classCol()<0 || from.classCol()>=from.dim() )
      throw whereOutOfRange(
	 "multiDS<T>::multiDS(const dataSet<T>&, int from.classCol())",
	 "from.classCol()",from.classCol(),0,from.dim()-1);
#endif
   // determine the number of distinct categories
   Matrix<int> singleton(1,1,from.classCol()); int garbage;
   dataSet<T> categories(from,singleton);
   categories.sort().checkMult(garbage);
   Matrix<T> catIndex(categories.matrix().ta());
   ds = new dataSet<T> [catIndex.n()]; label.resize(catIndex.n());
   for (int cI=0; cI<catIndex.n(); cI++)
   {
      Matrix<T> cm(from.distinct(),from.dim()-1);
      Matrix<int> multiple(from.distinct());
      Matrix<int> labelle(from.distinct());
      int cmPtr=0;
      for (int dPtr=0; dPtr<from.distinct(); dPtr++)
	 if ( from(dPtr,from.classCol())==catIndex(cI) )
	 {
	    if ( from.classCol() > 0 )
	       cm.s(cmPtr,1,0,from.classCol()) 
	       = from(dPtr).s(0,1,0,from.classCol());
	    if ( from.classCol() < from.dim()-1 )
	       cm.s(cmPtr,1,from.classCol()) 
	       = from(dPtr).s(0,1,from.classCol()+1);
	    multiple(cmPtr) = from.multiplicity(dPtr);
	    labelle(cmPtr) = from.label(dPtr);
	    cmPtr++;
	 }
      if ( cmPtr )
      {
	 ds[nbSets] = dataSet<T>(	cm.resize_(cmPtr,cm.n()),
				multiple.resize_(cmPtr),
				labelle.resize_(cmPtr));
	 label(nbSets) = catIndex(cI);
	 nbSets++;
      }
   }
   if ( nbSets < catIndex.n() ) // some categories were inexistant
   {
      label.resize_(nbSets);
      dataSet<T>* aux = ds;
      ds = new dataSet<T> [nbSets];
      for (int i=0; i<nbSets; i++) ds[i] = aux[i];
      delete [] aux;
   }
   for (int set=0; set<nbSets; set++)
      ds[set].monotone.resize(from.monotone);
   return *this;
}