/*
data structure of class cutPtsSet:
  private:
  Matrix<T>		origin;
  Matrix<float>	weight;
  Matrix<int>		var;
  Matrix<int>		sigmaCP;
  Matrix<float>	span;
  - *this is a row vector of nbCP() T, value of each cutPtsSet ***times 2***.
    the cutPtsSets are organized in an increasing order of the variables, and
    for each variable, in an increasing size.
  - weight is a row vector of nbCP() float, weight of each cutPtsSet.
  - var is a row vector of nbCP() int. var(cpIndex) is the index of the
    attribute associate to the cpIndex-th cutPtsSet.
  - sigmaCP is a row vector of nbVars() int.
    sigmaCP(var) := 	nbCP(var) + sigmaCP(var-1)	if var >  0
    nbCP(var)						if var == 0
  - span is a row vector of nbVars() float, span(var) is propotional to
    the maximal value along var, minus the minimal value along var.
  - origin is a row vector of nbCP() T, with the original value of each
    cutPtsSet. This is used only in case of normalization (by a call to
    normalize), otherwise this is an empty vector.
*/


#include "cutPtsSet.h"


#include <math.h>
#include <iostream.h>


#include "../Matrices/spareM.h"
#include "../Matrices/setCovering.h"
#include "main_def.h"



#define _CHECK_INDEX



// if T is char, short or int
#if 1
  #define Plus1(x) 	((x)+1)
  #define Div2(x) 	((x)>>1)
  #define Tim2(x) 	((x)<<1)
// else if T is float or double
#else
  #define Plus1(x) 	(x)
  #define Div2(x)	((x)/2.0)
  #define Tim2(x)	((x)*2.0)
#endif

								   


//////////////////////////// class cutPtsSet //////////////////////////////////

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

tcT cutPtsSet<T>::cutPtsSet (int nbCPts, int nbVars) : 
   Matrix<T>(nbCPts), origin(0),
   weight(nbCPts), var(nbCPts), sigmaCP(nbVars), span(nbVars)
{ }


tcT cutPtsSet<T>& cutPtsSet<T>::resize(int nbCPts, int nbVars)
{	Matrix<T>::resize(nbCPts); weight.resize(nbCPts); 
var.resize(nbCPts); sigmaCP.resize(nbVars); span.resize(nbVars);
if ( origin.n() ) origin.resize(nbCPts);
return *this;
}


tcT /*inline*/ cutPtsSet<T>::cutPtsSet () :
   Matrix<T>(0), origin(0), weight(0), var(0), sigmaCP(0), span(0)
{ }


tcT /*inline*/ cutPtsSet<T>::~cutPtsSet () { }


tcT cutPtsSet<T>::cutPtsSet (multiDS<T>& d, 
			     float& confidence, 
			     const char method/**=0**/, 
			     const char weighting/**=0**/) :
   Matrix<T>(0), origin(0), weight(0), var(0), sigmaCP(0), span(0)
{
   if ( method==0 )
      // insert a cut point at every change of class along the axis
   {
      resize(d.ds[0].dim()*d.distinct(), d.ds[0].dim());
      weight=0.0;
      Matrix<T> valcat(d.distinct(),2);
      int nbTh=0;
      for (int attr=0; attr<d.ds[0].dim(); attr++)
      {
	 sigmaCP(attr) = ( attr ? sigmaCP(attr-1) : 0);
	
	 int cl,clPtr;
	 for (cl=0, clPtr=0; cl<d.nbCats(); clPtr+=d.ds[cl].distinct(), cl++)
	 {
	    valcat.s(clPtr,d.ds[cl].distinct(),0,1) = d.ds[cl][attr];
	    valcat.s(clPtr,d.ds[cl].distinct(),1,1) = (T)(cl);
	 }
	 valcat.sort();  
	 int total_rows=valcat.m();
	 for(int first_d_k=0; 
	     first_d_k<total_rows && known(valcat(first_d_k,0)); 
	     first_d_k++);
	 if(first_d_k < total_rows)
	 {
	    int last_d_k;
	    for(last_d_k=first_d_k; 
		last_d_k+1<total_rows && unknown(valcat(last_d_k+1,0)); 
		last_d_k++);
	    valcat.s(first_d_k,total_rows-last_d_k-1) = 
	       valcat.s(last_d_k+1,total_rows-last_d_k-1);
	    total_rows -= last_d_k - first_d_k + 1;
	 }
	 span(attr) = float (valcat(total_rows-1,0) - valcat(0,0));
	 // Each row of valcat is a pair (value,class) in increasing order
	 // of the values. Many pairs can be identical, or with the same
	 // value but a different class.
	 // A cutPtsSet of value VAL is be introduced iff there is in valcat
	 // a row (VAL,C) and a row (VAL1,C1) such that C!=C1 and VAL1 is the
	 // largest possible value in valcat(.,0), strictly smaller than VAL.
	 boolean mixGroup=false; 
	 // true iff several rows in valcat, with current value VAL and 
	 // different classes, has been already identified
	 T prevThresh=(T)(valcat(0,0)); // value of the previous cutPtsSet
	 T val1; // VAL1 in above comment
	 boolean newCutPtsSet;
	 for (int i=1; i<total_rows; i++)
	 {
	    newCutPtsSet=false;
	    if ( valcat(i-1,0) != valcat(i,0) )
	    {
	       val1 = valcat(i-1,0);
	       newCutPtsSet = ( valcat(i-1,1)!=valcat(i,1) || mixGroup );
	       mixGroup=false;
	    }
	    else if ( valcat(i-1,1)!=valcat(i,1) ) // && vc(i-1,0)==vc(i,0)
	    {
	       newCutPtsSet = ( prevThresh != valcat(i,0) );
	       mixGroup=true;
	    }
	    if ( newCutPtsSet )
	    {
	       Array<T>::operator()(nbTh) = val1 + valcat(i,0);
	       if ( weighting==2 )
		  weight(nbTh) = float(valcat(i,0)-val1)/(2.0*span(attr));
	       var(nbTh) = attr;
	       sigmaCP(attr)++; 
	       nbTh++;
	       prevThresh=valcat(i,0);
	    }
	 }
      }
      resize_(nbTh); weight.resize_(nbTh); var.resize_(nbTh);
   }
   else if ( method==1 )
      // insert a cut point for every pairs of points of different classes
   {
      resize(d.ds[0].dim()*d.nbPairs(), d.ds[0].dim());
      Matrix<T> my_cut_pts(d.nbPairs());
      int nbTh=0;
      for (int attr=0; attr<d.ds[0].dim(); attr++)
      {
	 sigmaCP(attr) = ( attr ? sigmaCP(attr-1) : 0);
	 // compute the set of distinct values along attribute ATTR, for each
	 // class separatly
	 Matrix<int> singleton(1,1,attr);
	 multiDS<T> proj(d,singleton);
	 proj.sort().checkMult();
	 // introduce one cutPtsSet for each pairs of points from two
	 // different classes
	 int ptr=0; T ll,rr;
	 for (int clL=0; clL<proj.nbCats()-1; clL++)
	    for (int clR=clL+1; clR<proj.nbCats(); clR++)
	       for (int indL=0; indL<proj.ds[clL].distinct(); indL++)
		  for (int indR=0; indR<proj.ds[clR].distinct(); indR++)
		     if ( 
			(ll=proj.ds[clL](indL,0))!=(rr=proj.ds[clR](indR,0)) && 
			known(ll) && known(rr) &&
			  ( ll<rr ? d.monotone(attr)!=negativeAttr
			          : d.monotone(attr)!=positiveAttr ) 
			)
			my_cut_pts(ptr++) = ll+rr;
	  // sort and keep one cutPtsSet for each different value
	  if ( ptr )
	  {
	     my_cut_pts.s(0,1,0,ptr).sort();
	     T mini,maxi;
	     dataSet<T> merged;
	     merged = proj.merge(false);
	     merged.sort();
	     for(int i=0; i<merged.distinct() && 
		    unknown(merged(i,0)); i++);
	     if (i==merged.distinct()) span(attr)=1; 
	     else 
	     {
		mini=merged(i,0);
		for(i=merged.distinct()-1;
		    i>=0 && unknown(merged(i,0)); i--);
		maxi=merged(i,0);
		span(attr) = maxi - mini;
	     }
	     Array<T>::operator()(nbTh) = my_cut_pts(0);
	     var(nbTh) = attr; sigmaCP(attr)++; nbTh++;
	     for (int pth=1; pth<ptr; pth++)
		if ( my_cut_pts(pth)!=my_cut_pts(pth-1) )
		{
		   Array<T>::operator()(nbTh) = my_cut_pts(pth);
		   var(nbTh) = attr; sigmaCP(attr)++; nbTh++;
		}	
	  }
	  else span(attr)=0;		
      }
      resize_(nbTh); weight.resize_(nbTh); var.resize_(nbTh);
   }

   // Check whether confidence is not too high
   Array<int> list(d.listOfPairs());	
   float minConfi=1.0;
   for (int pair=0; pair<list.m(); pair++)
   {
      float maxGap=0.0;
      for (int pcp=0; pcp<nbCP(); pcp++)
      {
	 int vari=attribute(pcp);
	 T val1 = minimum(	d.ds[list(pair,0)](list(pair,1),vari),
				d.ds[list(pair,2)](list(pair,3),vari));
	 T val2 = maximum(	d.ds[list(pair,0)](list(pair,1),vari),
				d.ds[list(pair,2)](list(pair,3),vari));
	 if ( between(val1,val2,pcp) &&
	      ( val1<val2 ? d.monotone(vari)!=negativeAttr
		          : d.monotone(vari)!=positiveAttr ) && 
	      span(vari)>0.0 
	    )
	    maximize(maxGap,
		     float ((minimum(operator()(pcp)-Tim2(val1),
				    Tim2(val2)-operator()(pcp))) /
		      (2.0 * span(vari))));
      }
      minimize(minConfi,maxGap);
   }
   minConfi-=1E-5;
   if ( confidence > minConfi )
      flog<< endl << "The required confidence of " << setprecision(3) 
	 << setw(6) << confidence << " has been reset to " << setw(6)
	    << (confidence=minConfi) << endl << flush;
	

   // compute the weights
   if ( weighting==1 ) // correlation with output
   {
      for (int pCP=0; pCP<nbCP(); pCP++)
      {
	 int attr=attribute(pCP);
	 Matrix<double> cardi(2,d.nbCats(),0.0);
	 for (int set=0; set<d.nbCats(); set++)
	 {
	    for (int ind=0; ind<d.ds[set].distinct(); ind++)
	       switch ( compare(d.ds[set](ind,attr),pCP,confidence) )
	       {
		case +1: { 	cardi(1,set)+=d.ds[set].multiplicity(ind); 
		break; 
		}
		case  0:	cardi(0,set)+=d.ds[set].multiplicity(ind);
	       }
	 }
	 if ( cardi.sum().min()(0) = 0.0 )
	    // this can never happen when confidence==0. However, when
	    // confidence is >0, it might happen in some rare occasions
	    weight(pCP) = -1E10;
	 else
	 {
	    cardi /= cardi.sum().t();
	    double entropyPos=0.0; double entropyNeg=0.0;
	    for (set=0; set<d.nbCats(); set++)
	    {
	       entropyPos += ( cardi(1,set)>0.0 ? 
			       cardi(1,set)*log(cardi(1,set)) : 0.0 );
	       entropyNeg += ( cardi(0,set)>0.0 ? 
			       cardi(0,set)*log(cardi(0,set)) : 0.0 );
	    }
	    weight(pCP) = float(maximum(entropyPos,entropyNeg));
	 }
      }
   }
   else if ( weighting==2 && method )
   {
      weight=1E10;
      for (int attr=0; attr<d.ds[0].dim(); attr++)
      {
	 Matrix<int> singleton(1,1,attr);
	 multiDS<T> proj(d,singleton);
	 proj.sort().checkMult();
	 for (int clL=0; clL<proj.nbCats()-1; clL++)
	    for (int clR=clL+1; clR<proj.nbCats(); clR++)
	       for (int indL=0; indL<proj.ds[clL].distinct(); indL++)
		  for (int indR=0; indR<proj.ds[clR].distinct(); indR++)
		  {
		     int val1 = proj.ds[clL](indL,0);
		     int val2 = proj.ds[clR](indR,0);
		     if ( val1 != val2 && known(val1) && known(val2))
		     {
			int first=firstCPG( minimum(val1,val2),
					    attr,confidence);
			int last = lastCPL( maximum(val1,val2),
					    attr,confidence);
			if ( first>=0 )
			   for (int pth=first; pth<=last; pth++)
			      minimize(weight(pth),float(minimum(
				 absolute(Tim2(val1)-operator()(pth)),
				 absolute(Tim2(val2)-operator()(pth))))*
				       float(proj.ds[clL].multiplicity(indL))*
				       float(proj.ds[clR].multiplicity(indR)));
		     }
		  }
	 for (int pth=firstIndex(attr); 
	      pth<firstIndex(attr)+nbCP(attr); pth++)
	    if ( weight(pth)<1E10 )
	       weight(pth)/=(2.0*span(attr));
	    else
	       weight(pth)=0.0;
      }
   }
   else if ( weighting==3 ) // global separability
   {   for (int attr=0; attr<d.ds[0].dim(); attr++)
   {
      Matrix<int> singleton(1,1,attr);
      multiDS<T> proj(d,singleton);
      proj.sort().checkMult();
      for (int clL=0; clL<proj.nbCats()-1; clL++)
	 for (int clR=clL+1; clR<proj.nbCats(); clR++)
	    for (int indL=0; indL<proj.ds[clL].distinct(); indL++)
	       for (int indR=0; indR<proj.ds[clR].distinct(); indR++)
	       {
		  int val1 = proj.ds[clL](indL,0);
		  int val2 = proj.ds[clR](indR,0);
		  if ( val1 != val2 && known(val1) && known(val2))
		  {
		     int first=firstCPG( minimum(val1,val2),
					 attr,confidence);
		     int last = lastCPL( maximum(val1,val2),
					 attr,confidence);
		     if ( first>=0 )
			for (int pth=first; pth<=last; pth++)
			   weight(pth) += float(minimum(
			      absolute(Tim2(val1)-operator()(pth)),
			      absolute(Tim2(val2)-operator()(pth))))*
			      float(proj.ds[clL].multiplicity(indL))*
			      float(proj.ds[clR].multiplicity(indR));
		  }
	       }
      weight.s(0,1,firstIndex(attr),nbCP(attr))/=(2.0*span(attr));
   }
   }
}


tcT cutPtsSet<T>::cutPtsSet (cutPtsSet<T>& from, const Matrix<int>& index) :
   Matrix<T>(index.n()), origin(0),
   weight(index.n()), var(index.n()), 
   sigmaCP(1,from.nbVars(),0), span(from.span)
{
   if ( from.origin.n() ) origin.resize(index.n());
   Matrix<int> ind(index); ind.sort();
   for (int th=0; th<ind.n(); th++)
   {
      Array<T>::operator()(th) = from(ind(th));
      if ( origin.n() ) origin(th) = from.origin(ind(th));
      weight(th) = from.weight(ind(th));
      var(th) = from.var(ind(th));
      if ( th && var(th)!=var(th-1) )
	 for (int j=var(th-1)+1; j<=var(th); j++)
	    sigmaCP(j) = sigmaCP(j-1);
      sigmaCP(var(th))++;
   }
   for (int j=var(th-1)+1; j<sigmaCP.n(); j++)
      sigmaCP(j) = sigmaCP(j-1);
}



///////////////
// Selectors //
///////////////

tcT /*inline*/ int cutPtsSet<T>::nbVars () const 
{ return sigmaCP.n(); }


tcT /*inline*/ int cutPtsSet<T>::nbCP () const 
{ return sigmaCP(nbVars()-1); }


tcT /*inline*/ int cutPtsSet<T>::nbCP (int vari) const
{
#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	"int cutPtsSet<T>::nbCP(int vari) const",
				"vari", vari, 0, nbVars()-1);
#endif
   if ( !vari ) return sigmaCP(0);
   else return sigmaCP(vari)-sigmaCP(vari-1);
}


tcT /*inline*/ int cutPtsSet<T>::firstIndex (int vari) const
{
#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	"int cutPtsSet<T>::firstIndex(int vari) const",
				"vari", vari, 0, nbVars()-1);
#endif
   return (vari ? sigmaCP(vari-1) : 0); 
}


tcT /*inline*/ int cutPtsSet<T>::lastIndex  (int vari) const 
{
#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	"int cutPtsSet<T>::lastIndex(int vari) const",
				"vari", vari, 0, nbVars()-1);
#endif
   return sigmaCP(vari)-1; 
}



tcT /**inline**/ int cutPtsSet<T>::attribute (int cpIndex) const
{
#ifdef _CHECK_INDEX
   if ( cpIndex<0 || cpIndex>=nbCP() )
      throw whereOutOfRange(	"int cutPtsSet<T>::attribute(int thInd) const",
				"thInd", cpIndex, 0, nbCP()-1);
#endif
   return var(cpIndex);
}


tcT /*inline*/ int cutPtsSet<T>::index (int cpIndex) const
{
#ifdef _CHECK_INDEX
   if ( cpIndex<0 || cpIndex>=nbCP() )
      throw whereOutOfRange(	"int cutPtsSet<T>::index(int cpIndex) const",
				"cpIndex", cpIndex, 0, nbCP()-1);
#endif
   return (var(cpIndex) ? cpIndex - sigmaCP(var(cpIndex)-1) : cpIndex); 
}


tcT /*inline*/ T cutPtsSet<T>::operator () (int cpIndex) const
{

#ifdef _CHECK_INDEX
   if ( cpIndex<0 || cpIndex>=nbCP() )
      throw whereOutOfRange(	"T cutPtsSet<T>::operator()(int cpIndex) const",
				"cpIndex", cpIndex, 0, nbCP()-1);
#endif
   return Array<T>::operator()(cpIndex); 
}


tcT /**inline**/ boolean cutPtsSet<T>::above (T val, int cpIndex) const
{
   if(unknown(val)) return(false);
#if _val_GEQ_cp
   return ( Tim2(val) >= operator()(cpIndex) ); 
#else
   return ( Tim2(val) >  operator()(cpIndex) );
#endif
}


tcT int cutPtsSet<T>::compare (T val, int cpIndex, 
			     float confidence/**=0.0**/) const
{
   if(unknown(val)) return(-1);	
   confidence *= span(attribute(cpIndex));
#if _val_GEQ_cp
   if ( float(val) >= float(operator()(cpIndex))/2.0 + confidence )
      return +1;
   else if ( float(val) < float(operator()(cpIndex))/2.0 - confidence )
      return  0;
   else
      return -1;
#else
   if ( float(val) > float(operator()(cpIndex))/2.0 + confidence )
      return +1;
   else if ( float(val) <= float(operator()(cpIndex))/2.0 - confidence )
      return  0;
   else
      return -1;
#endif
}


tcT	boolean cutPtsSet<T>::between (T val1, T val2, int cpIndex, 
				     float confidence/**=0.0**/) const
{
   if ( val1==val2 || unknown(val1) || unknown(val2)) return false;
   confidence *= span(attribute(cpIndex));
#if _val_GEQ_cp
   if ( val1 < val2 ) 
      return float(val1)+confidence <  float(operator()(cpIndex))/2.0
      && float(val2)-confidence >= float(operator()(cpIndex))/2.0;
   else
      return float(val2)+confidence <  float(operator()(cpIndex))/2.0
      && float(val1)-confidence >= float(operator()(cpIndex))/2.0;
#else
   if ( val1 < val2 ) 
      return float(val1)+confidence <= float(operator()(cpIndex))/2.0
      && float(val2)-confidence >  float(operator()(cpIndex))/2.0;
   else
      return float(val2)+confidence <= float(operator()(cpIndex))/2.0
      && float(val1)-confidence >  float(operator()(cpIndex))/2.0;
#endif
}


tcT /*inline*/ T cutPtsSet<T>::operator () (int cpInV, int vari) const
{

#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::operator()(int cpInV, int vari) const",
	 "vari", vari, 0, nbVars()-1);
   if ( cpInV<0 || cpInV>=nbCP(vari) )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::operator()(int cpInV, int vari) const",
	 "cpInV", cpInV, 0, nbCP(vari)-1);
#endif
   return Array<T>::operator()(firstIndex(vari)+cpInV); 
}


tcT /*inline*/ boolean cutPtsSet<T>::above (T val, int cpInV, int vari) const
{
   if(unknown(val)) return(false);
#if _val_GEQ_cp
   return ( Tim2(val)>=operator()(cpInV,vari) ); 
#else
   return ( Tim2(val)> operator()(cpInV,vari) );
#endif
}


tcT int cutPtsSet<T>::compare (T val, int cpInV, int vari, 
			     float confidence/**=0.0**/) const
{
   if(unknown(val)) return(-1);		
   confidence *= span(vari);
#if _val_GEQ_cp
   if ( float(val) >= float(operator()(cpInV,vari))/2.0 + confidence )
      return +1;
   else if ( float(val) < float(operator()(cpInV,vari))/2.0 - confidence )
      return -1;
   else
      return 0;
#else
   if ( float(val) > float(operator()(cpInV,vari))/2.0 + confidence )
      return +1;
   else if ( float(val) <= float(operator()(cpInV,vari))/2.0 - confidence )
      return -1;
   else
      return 0;
#endif
}


tcT	boolean cutPtsSet<T>::between (T val1, T val2, int cpInV, int v,
				     float confidence/**=0.0**/) const
{
   if ( val1==val2 || unknown(val1) || unknown(val2) ) return false;
#if _val_GEQ_cp
   if ( val1 < val2 ) 
      return float(val1)+confidence <  float(operator()(cpInV,v))/2.0
      && float(val2)-confidence >= float(operator()(cpInV,v))/2.0;
   else
      return float(val2)+confidence <  float(operator()(cpInV,v))/2.0
      && float(val1)-confidence >= float(operator()(cpInV,v))/2.0;
#else
   if ( val1 < val2 ) 
      return float(val1)+confidence <= float(operator()(cpInV,v))/2.0
      && float(val2)-confidence >  float(operator()(cpInV,v))/2.0;
   else
      return float(val2)+confidence <= float(operator()(cpInV,v))/2.0
      && float(val1)-confidence >  float(operator()(cpInV,v))/2.0;
#endif
}


tcT /*inline*/ T cutPtsSet<T>::value (int cpInV, int vari) const
{

#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::operator()(int cpInV, int vari) const",
	 "vari", vari, 0, nbVars()-1);
   if ( cpInV<0 || cpInV>=nbCP(vari) )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::operator()(int cpInV, int vari) const",
	 "cpInV", cpInV, 0, nbCP(vari)-1);
#endif
#if _val_GEQ_cp
   return Div2(Plus1(Array<T>::operator()(firstIndex(vari)+cpInV))); 
#else
   return Div2(Array<T>::operator()(firstIndex(vari)+cpInV)); 
#endif
}


tcT /*inline*/ float cutPtsSet<T>::original (int cpInV, int vari) const
{

#ifdef _CHECK_INDEX
   if ( vari<0 || vari>=nbVars() )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::original(int cpInV, int vari) const",
	 "vari", vari, 0, nbVars()-1);
   if ( cpInV<0 || cpInV>=nbCP(vari) )
      throw whereOutOfRange(	
	 "T cutPtsSet<T>::original(int cpInV, int vari) const",
	 "cpInV", cpInV, 0, nbCP(vari)-1);
#endif
   return ( origin.n() ? float(origin(firstIndex(vari)+cpInV))/2.0
	    : float(Array<T>::operator()(firstIndex(vari)+cpInV))/2.0); 
}


tcT int cutPtsSet<T>::firstCPG (T val, int vari, float confidence/**=0.0**/) const
{
   if(unknown(val)) return -1;	
   int last = lastIndex(vari);
   for (int index=firstIndex(vari); index<=last; index++)
      if ( compare(val,index,confidence)==0 ) return index;
   return -1;
}

tcT int cutPtsSet<T>::lastCPL (T val, int vari, float confidence/**=0.0**/) const
{
   if(unknown(val)) return -1;	
   int first = firstIndex(vari);
   for (int index=lastIndex(vari); index>= first; index--)
      if ( compare(val,index,confidence)==1 ) return index;
   return -1;
}



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

tcT cutPtsSet<T>& cutPtsSet<T>::normalize (multiDS<T>& d)
{
   if ( d.normalized ) return *this;
   for (int set=0; set<d.nbCats(); set++)
   {
      Matrix<T> nd(d.ds[set].distinct(),d.ds[set].dim());
      nd = d.ds[set].matrix();
      for (int vari=0; vari<nd.n(); vari++)
	 for (int row=0; row<nd.m(); row++)
	    if ( known(nd(row,vari)) ) 
	       for (int tI=0; tI<nbCP(vari); )
		  if ( !above(nd(row,vari),tI,vari) || ++tI>=nbCP(vari) )
		  {
		     nd(row,vari) = (T)(tI);
		     break; 
		  }
      d.ds[set] = dataSet<T>(nd,d.ds[set]);
   }
   if ( !origin.n() ) origin.resize(*this);
   for (int row=0; row<nbVars(); row++)
      for (int tI=0; tI<nbCP(row); tI++)
	 Array<T>::operator()(firstIndex(row)+tI) = (T)(2*tI+1);
   d.normalized=true;
   return *this;
}


tcT cutPtsSet<T>& cutPtsSet<T>::reduce (multiDS<T>& d,
					int method, 
					float& separability, 
					float confidence/**=0.0**/)
{
   // Determine the set of pairs of points from two different classes.
   Array<int> list(d.listOfPairs());	

   // The indices of the final cut points will be ssind(0),...,ssind(ssSize-1).
   Matrix<int> ssind(0); int ssSize=0;


   if ( method<0 )
      cout<< endl << "There are X cut points of weight>=Y, of min sep Z\n";


   if ( method<=1 )
   {
      ssind.resize(1,nbCP()); 
      for (int i=0; i<nbCP(); i++) ssind(i)=i;
      weight.sort(ssind);
      Matrix<float> separabil(1,list.m()); separabil=0.0; float sepMin;
      for (int curW=nbCP()-1; curW>=0; curW--)
      {
	 int vari=attribute(ssind(curW));
	 for (int pair=0; pair<list.m(); pair++)
	 {
	    T val1=d.ds[list(pair,0)](list(pair,1),vari);
	    T val2=d.ds[list(pair,2)](list(pair,3),vari);		
	    if ( between(val1,val2,ssind(curW),confidence) &&
		 ( val1 < val2 ? d.monotone(vari)!=negativeAttr
		               : d.monotone(vari)!=positiveAttr )
	       )
	       if ( method%2 )
		  if ( val1 < val2 )
		     separabil(pair) += float (
			minimum(Tim2(val2)-operator()(ssind(curW)),
				operator()(ssind(curW))-Tim2(val1)))
		     / (2.0 * span(vari));
		  else
		     separabil(pair) += float (
			minimum(Tim2(val1)-operator()(ssind(curW)),
				operator()(ssind(curW))-Tim2(val2)))
		     / (2.0 * span(vari));
	       else
		  separabil(pair) += 1.0;
	 }
	 ssSize++;
	 sepMin = separabil.min()(0);
	 if ( method<0 )
	 {
	    if ( !curW || weight(ssind(curW))!=weight(ssind(curW-1)) )
	       cout<< setw(11) << ssSize << setw(24) << setprecision(5)
		  << weight(ssind(curW)) << setw(14) << sepMin << endl;
	    if ( !curW ) return *this;
	 }
	 else if ( sepMin >= separability )
	 {
	    while ( curW && weight(ssind(curW))==weight(ssind(--curW)) )
	       ssSize++;
	    break;
	 }
      }
      if ( ssSize < nbCP() ) 
	 ssind.s(0,1,0,ssSize) = ssind.s(0,1,nbCP()-ssSize,ssSize);
   }
	

   if ( method==2 )
   {
      binMatrix constraint(list.m(), nbCP(), false);
      for (int row=0; row<list.m(); row++)
      {
	 int first, last;
	 int val1, val2;
	 int totSep=0;
	 for (int vari=0; vari<nbVars(); vari++)
	 {
	    val1 = d.ds[list(row,0)](list(row,1),vari);
	    val2 = d.ds[list(row,2)](list(row,3),vari);
	    if ( val1 != val2 && known(val1) && known(val2) &&
	         ( val1<val2 ? d.monotone(vari)!=negativeAttr
		             : d.monotone(vari)!=positiveAttr )
	       )
	    {
	       first = firstCPG(minimum(val1,val2),vari,confidence);
	       last  = lastCPL (maximum(val1,val2),vari,confidence);
	       if ( first>=0 && last>=0 && last>=first )
	       {
		  constraint.s(row,1,first,last-first+1)=true;
		  totSep+=last-first+1;
	       }
	    }
	 }

	 if ( totSep < separability )
	    flog<< endl << "### Pair (" 
		<< d.ds[list(row,0)].label(list(row,1)) << ","
		<< d.ds[list(row,2)].label(list(row,3)) 
		<< ") has a total separability of" 
		<< setw(3) << totSep << ", strictly less than " 
		<< setw(3) << separability << endl << flush;
	 else if ( !totSep )
	    throw whereWhich(
	       "cutPtsSet<T>::cutPtsSet (const multiDS<T>& d, float& SEP, int)",
	       "Unseparable pair of points");
      }

      setCover sc(constraint); 
      sc.greedy(flog, 1,int(separability),false); 
      if ( int(separability)==1 ) sc.reduce();
      ssind.resize(sc.solution()); ssSize=ssind.n();
   }


   if ( method==3 )
   {
      // flog << "\n Hello, I am entering the greedy GSCP\n" << flush;
      // Construct a matrix SCP with one column per pair of points from 
      // differents classes, and one row per current cut points, such that 
      // SCP(k,(i,j)) gives a non-negative value expressing the quality of the
      // separation of the points i and j using cutPtsSet k: SCP(k,(i,j)) = 0
      // if k does not separate i and j, and than SCP(k,(i,j)) is proportional
      // to the minimum gap between i and k and j and k.
      // In this implementation, SCP(k,(i,j)) is in {0,...,60}. 60 corresponds
      // to a cutPtsSet of value 0.5 for a pair of boolean points (0,1).

      // Row 0 and columns 0-1 of SCP will be used later.
      const int hOff=2; const int vOff=1; // col and row offsets

      int bound = int(separability * 120.0);
      // Thus separability==0.5 implies bound==60.
      if ( bound > 240 ) // this will lead to results of sums above 255
	 throw whereWhich(
	    "cutPtsSet<T>::cutPtsSet (const multiDS<T>& d, float& SEP, int)",
	    "SEP can not be bigger than 2");

      // flog << "There are " << list.m() << " pairs and " << nbCP() << " cut points\n";
      // flog << "I try to declare SCP\n" << flush;

      // Count the number of non zero coefficient of the matrix
      int nbNZ=0;
      int pair;
      for (pair=0; pair<list.m(); pair++)
      {
	 for (int th=0; th<nbCP(); th++)
	 {
	    int vari=attribute(th);
	    T val1=d.ds[list(pair,0)](list(pair,1),vari);
	    T val2=d.ds[list(pair,2)](list(pair,3),vari);		
	    if ( between(val1,val2,th,confidence) &&
		 ( val1 < val2 ? d.monotone(vari)!=negativeAttr
		               : d.monotone(vari)!=positiveAttr )
	       )
	       nbNZ++;
	 }
      }

      // declare the constraint matrix
      spareM<unsigned char> SCP(nbCP(),nbNZ); 

      // setup the constraint matrix and compute the total sum for each pair
      Matrix<int> sumTot(list.m());
      unsigned char aux;
      for (int th=0; th<nbCP(); th++)
      {
	 int vari=attribute(th);
	 for (pair=0; pair<list.m(); pair++)
	 {
	    T val1=d.ds[list(pair,0)](list(pair,1),vari);
	    T val2=d.ds[list(pair,2)](list(pair,3),vari);		
	    if ( between(val1,val2,th,confidence) &&
		 ( val1 < val2 ? d.monotone(vari)!=negativeAttr
		               : d.monotone(vari)!=positiveAttr )
	       )
	    {
	       if ( val1 < val2 )
		  aux = ceil(60.0*float(minimum(
		     Tim2(val2)-operator()(th),
		     operator()(th)-Tim2(val1)))/span(vari));
	       else
		  aux = ceil(60.0*float(minimum(
		     Tim2(val1)-operator()(th),
		     operator()(th)-Tim2(val2)))/span(vari));
	       SCP.newCoef(th,pair,aux);
	       sumTot(pair) += aux;
	    }
	 }
      }

      // Find out the critical pairs
      Matrix<int> criticalPairs(list.m()); int nbcriticalPairs=0;
      for (pair=0; pair<list.m(); pair++)
	 if ( sumTot(pair) <= bound )
	 {
	    criticalPairs(nbcriticalPairs++)=pair;
	    if ( sumTot(pair) < bound )
	       flog<< endl << "### Pair (" 
		   << d.ds[list(pair,0)].label(list(pair,1)) << ","
		   << d.ds[list(pair,2)].label(list(pair,3)) 
		   << ") has a total separability of " 
		   << setw(8) << sumTot(pair)/120.0 << ", strictly less than " 
		   << setw(8) << bound/120.0 << endl << flush;
	 }

      // Solve the ILP problem:  min sum Xi, st SCP^t*X >= bound, Xi in {0,1}.
      // row 0 of SCP contains the sum due to the current solution
      // col 0 of SCP contains the index of the cutPtsSet
      ssind.resize(nbCP()); ssSize=0;
      int nbR=nbCP(); int nbC=list.m();
      int row,col; unsigned char v;
      int bestRow; int minToFill; int toFill;
      Matrix<int> auxsum(1,list.m(),0);
      Matrix<int> partsum(1,list.m(),0);
      Matrix<char> activeR(1,nbCP(),true);
      Matrix<char> activeC(1,list.m(),true);
      Matrix<float> lweight(weight);

      // if there are some pairs with separability<=required_separability
      // then we first select all the cut points that separate these pairs
      if ( nbcriticalPairs )
      {
	 // Introduce all the required rows in the solution, suppressed them
	 // from the table, and mark all the columns sufficiently covered
	 for (int pairInd=0; pairInd<nbcriticalPairs; pairInd++)
	 {
	    activeC(criticalPairs(pairInd))=false; nbC--;
	    for (row=0; row<nbCP(); row++)
	       if ( activeR(row) && SCP(row,criticalPairs(pairInd)) )
	       {
		  activeR(row)=false; nbR--;
		  ssind(ssSize++) = row;
		  if ( SCP.first(row,col,v) )
		     do
			partsum(col) += v;
		     while ( SCP.next(row,col,v) );
	       }
	 }
	 // delete all the columns sufficiently covered
	 for (int cc=0; cc<list.m(); cc++)
	    if ( activeC(cc) && partsum(cc)>=bound ) 
	    {
	       activeC(cc)=false; nbC--;
	    }
      }

      // flog << "I am starting solving the problem\n" << flush;
      // main loop
      while (nbC && nbR)
      {
	 // find the best row
	 minToFill=nbC*bound; bestRow=0;
	 for (row=0; row<nbCP(); row++)
	    if ( activeR(row) )
	    {
	       auxsum = partsum;
	       if ( SCP.first(row,col,v) )
		  do
		     if ( activeC(col) ) auxsum(col) += v;
		  while ( SCP.next(row,col,v) );
	       
	       for (col=0, toFill=0; col<list.m(); col++)
		  if ( activeC(col) && auxsum(col) < bound ) 
		     toFill+=(bound-auxsum(col));
	       if ( (toFill<minToFill) 
		    || ( (toFill==minToFill) 
			 && ( (lweight(row)>lweight(bestRow))
			      || ((lweight(row)==lweight(bestRow)) && flipAtail)
			    )
		       )
		  ) 
	       {
		  minToFill=toFill; bestRow=row;
	       }
	    }
	 // add bestRow into solution, update the sum and suppress bestRow
	 ssind(ssSize++) = bestRow;
	 activeR(bestRow)=false; nbR--;
	 if ( SCP.first(bestRow,col,v) )
	    do
	       if ( activeC(col) ) 
		  if ( partsum(col)+v < bound )
		     partsum(col) += v;
		  else
		  { activeC(col)=false; nbC--; }
	    while ( SCP.next(bestRow,col,v) );

      }
   }
   // flog << "I did it!\n" << flush;


   // Decode the solution, and construct the final set of cut points.

   ssind.s(0,1,0,ssSize).sort();
   // Since ssind represents a subset and the indices a orders, one can
   // assume that ssind(i)>=i for all i=0,...,ssSize-1.
   sigmaCP = 0;
   for (int s=0; s<ssSize; s++)
   {
      Array<T>::operator()(s) = Array<T>::operator()(ssind(s));
      var(s) = var(ssind(s));
      weight(s) = weight(ssind(s));
      if ( origin.n() ) origin(s) = origin(ssind(s));
      if ( s && var(s)!=var(s-1) ) 
	 sigmaCP.s(0,1,var(s-1)+1,var(s)-var(s-1)) = sigmaCP(var(s-1));
      sigmaCP(var(s))++;
   }
   if ( var(s-1)+1 < sigmaCP.n() )
      sigmaCP.s(0,1,var(s-1)+1) = sigmaCP(var(s-1));
   resize_(ssSize); weight.resize_(ssSize); var.resize_(ssSize);
   if ( origin.n() ) origin.resize_(ssSize);
   return *this;
}


tcT ostream& operator << (ostream& s, const cutPtsSet<T>& myself)
{
   s.setf(ios::fixed | ios::showpoint);
   s << endl;
   int curCP=0;
   for (int v = 0; v < myself.nbVars(); v++)
      if ( myself.nbCP(v) )
      {
	 s 	<< "v" << setw(2) << v+1 << ":";
	 for (int ptr=0; ptr<myself.nbCP(v); ptr++) 
	 {
	    if ( ptr && !(ptr%5) ) s << endl << "    ";
	    s	<< setw(6) << ++curCP << ":"
		<< setw(6) << setprecision(1) << myself.original(ptr,v)
		   << setw(7) << setprecision(3) 
		      << myself.weight(myself.firstIndex(v)+ptr);
	 }
	 s << endl;
      }
   return s << "#total :" << myself.nbCP() << endl;
}