dataGen.cc

#include<iostream>
#include<cmath>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<fstream>
#include<algorithm>
#include<ext/algorithm> // for linux code
//#include<algorithm>     // for solaris/unix code
#include<utility>
#include<unistd.h>      // for getopt() function
using namespace std;
using namespace __gnu_cxx; // for linux code

enum Digit {MINUS_ONE=-1,ZERO,ONE,TWO,THREE,FOUR,FIVE,SIX,SEVEN,EIGHT,NINE,TEN};
enum DistMode {FULL_DIM=0,HYPER_RECT_UNIFORM,MVNORMAL,LINEAR_DEPENDENCE};
static const int RANGE = 1000;
static const int SAMPLE_SIZE = 1000000;
static const int NUM_ARGS = TWO;
//input args
int iNoOfClusters=5,iNoOfPts=1000,iNoOfDims=400,iRealDims=50,iMode=2,iNoOfDatasets=1;
double dProbNoise=0.05,dStdDev=0.5,dMaxClusterSzRatio=4.0,dProbOverlap = 0.5,dProbConstrain=0.5,x = 0.5,dBias = 0.05,dWidth=20;//0.25*RANGE;
char caOutFile[300];
bool bCentresOnly = false;

// see Probability Methods for Computer Science - Sheldon Ross (Chpt 9)
double genStdNormal() {   
  while(1) {
    double dY1 = -log(1.0 - (rand()%RANGE)/(double)RANGE);
    double dY2 = -log(1.0 - (rand()%RANGE)/(double)RANGE);
    double dY = (dY2-(dY1-1)*(dY1-1))/2.0;
    if(dY >= 0.0) return ((rand()%RANGE > RANGE/TWO) ? dY1 : -dY1);
  }
}

double estimateHashProb(double dStdDev,int iDiv=10) {
  int iNoOfHits = 0;
  double dDiv = 1.0/iDiv;
  for(int i=0;i<SAMPLE_SIZE;++i) {
    double dY = double(rand()%(int(dDiv*RANGE)))/RANGE;
    double dSample = genStdNormal();
    if(dSample <= (-dY+dDiv)/dStdDev && dSample >= -dY/dStdDev) 
      ++iNoOfHits;
  }
  return double(iNoOfHits)/SAMPLE_SIZE;
}

void createClusterSizes(int* iaCoord,double dBias,int iNoOfClusters,double v,double dProbNoise,ofstream& fOut,int RANGE=1000) {
  int iaRange[RANGE];
  for(int i = MINUS_ONE;++i < RANGE;) iaRange[i] = i;
  // generate P(i \in jth cluster) randomly with uniform bias dBias
  int iMinMinAlpha = int(max(dBias,(1.0-dProbNoise)/(ONE+v*(iNoOfClusters-ONE)))*RANGE);
  int iMaxMinAlpha = int(RANGE*(1.0-dProbNoise)/(v+iNoOfClusters-ONE));
  if(iMaxMinAlpha > iMinMinAlpha) {
    iaCoord[ZERO]=rand()%(iMaxMinAlpha-iMinMinAlpha)+iMinMinAlpha;
  } else iaCoord[ZERO]=iMinMinAlpha;
  int iRangeSize = int(RANGE*(1.0-dProbNoise)-iaCoord[ZERO]*(iNoOfClusters-TWO+v+ONE));
  if(iRangeSize<iNoOfClusters-TWO) {
    cerr << "ERROR: UNSUITABLE PARAMETERS: Decrease -b " << dBias << " or -a " << v << endl;
    exit(0);
  }
  random_sample_n(iaRange,iaRange+iRangeSize,iaCoord+ONE,iNoOfClusters-TWO);
  for(int j=ZERO;++j < iNoOfClusters-TWO;) iaCoord[j]=iaCoord[j+1];
  iaCoord[iNoOfClusters-TWO]=iRangeSize;
  sort(iaCoord+ONE,iaCoord+iNoOfClusters-ONE);
  for(int i=ZERO;++i<iNoOfClusters-ONE;)iaCoord[i] += ((i+1)*iaCoord[ZERO]);
  iaCoord[iNoOfClusters-ONE]=int((1.0-dProbNoise)*RANGE);
  int iTmp = ZERO;
  fOut << iNoOfClusters << ' ' << dProbNoise << ' ';
  for(int i=0;i < iNoOfClusters;++i) {
    fOut << double(iaCoord[i]-iTmp)/RANGE << ' ';
    iTmp=iaCoord[i];
  }
  fOut << endl;
}

void createClusterParams(int* iaMean,double *daStdDev,int iNoOfDims,int iNoOfClusters,DistMode distMode,int iDims,double dProbConstrain=0.5,double dProbOverlap=0.5) {
  int iaRange[RANGE];
  double r=2;
  int s=2;
  for(int iCntr = MINUS_ONE;++iCntr < iNoOfClusters;) {
    for(int d = MINUS_ONE;++d < iNoOfDims;) {
/*      int iScaleFactor = (rand()%s) + ONE;
      if(distMode==MVNORMAL || distMode==HYPER_RECT_UNIFORM)
      daStdDev[iNoOfDims*iCntr+d]=iScaleFactor*iScaleFactor*r*r;
      daStdDev[iNoOfDims*iCntr+d]=double(10+rand()%20);*/
      daStdDev[iNoOfDims*iCntr+d]=dWidth;
    }
  }
  for(int i = MINUS_ONE;++i < RANGE;) iaRange[i] = i;
  double x = 0.8;
  int iaCoord[iNoOfClusters];

  // generate the random full-dimensional centres
  switch(distMode) {
    case FULL_DIM:
      for(int d = MINUS_ONE;++d < iNoOfDims;) {
	// choose iNoOfClusters random points in range of dimension as centres
        random_sample_n(iaRange,iaRange+RANGE,iaCoord,iNoOfClusters);
        random_shuffle(iaCoord,iaCoord+iNoOfClusters);
        for(int iCntr = MINUS_ONE;++iCntr < iNoOfClusters;) {
          iaMean[iNoOfDims * iCntr + d] = iaCoord[iCntr];
        }
      }
      break;

    case HYPER_RECT_UNIFORM:
      for(int iCntr = MINUS_ONE;++iCntr < iNoOfClusters;) {
	double c = dProbConstrain*((iNoOfDatasets > 1)?(x+ (rand()%int((1.0-x)*1000))/500.0):1);
	double o = dProbOverlap*((iNoOfDatasets > 1)?(x+ (rand()%int((1.0-x)*1000))/500.0):1);
        double dO = ((iCntr > ZERO) ? (1.0-o)/(1.0-c) : 1.0); //~
        for(int d = MINUS_ONE;++d < iNoOfDims;) {
          if(iCntr > ZERO && iaMean[iNoOfDims*(iCntr-1)+d]!=MINUS_ONE) {
            if(rand()%100 < o*100)
              iaMean[iNoOfDims*iCntr+d]=iaMean[iNoOfDims*(iCntr-1)+d]+2*int(daStdDev[iNoOfDims*(iCntr-1)+d]*(((rand()%2)==0)?-1:1));//2 \sigma from prev mean
          } else {
            if(rand()%100 < c*dO*100)
              iaMean[iNoOfDims*iCntr+d]=rand()%RANGE;
          }
        }
      }
      break;

    case MVNORMAL:
      for(int iCntr = MINUS_ONE;++iCntr < iNoOfClusters;) {
	double c = dProbConstrain*((iNoOfDatasets > 1)?(x+ (rand()%int((1.0-x)*1000))/500.0):1);
	double o = dProbOverlap*((iNoOfDatasets > 1)?(x+ (rand()%int((1.0-x)*1000))/500.0):1);
        double dO = ((iCntr > ZERO) ? (1.0-o)/(1.0-c) : 1.0);
        for(int d = MINUS_ONE;++d < iNoOfDims;) {
          if(iCntr > ZERO && iaMean[iNoOfDims*(iCntr-1)+d]!=MINUS_ONE) {
            if(rand()%100 < o*100)
              iaMean[iNoOfDims*iCntr+d]=iaMean[iNoOfDims*(iCntr-1)+d]+2*int(daStdDev[iNoOfDims*(iCntr-1)+d]*(((rand()%2)==0)?-1:1));//2 \sigma from prev mean
//              iaMean[iNoOfDims*iCntr+d]=iaMean[iNoOfDims*(iCntr-1)+d]+(rand()%(2*int(daStdDev[iNoOfDims*(iCntr-1)+d])))*((rand()%2)?-1:1);
          } else {
            if(rand()%100 < c*dO*100)
              iaMean[iNoOfDims*iCntr+d]=rand()%RANGE;
          }
        }
      }
      break;

    default:
      cerr << "ERROR: UNDEFINED DATA GENERATION MODE " << distMode << endl;
      exit(0);
  }
#ifdef OLD_DEBUG
  for(int i=MINUS_ONE;++i < iNoOfClusters;) {
    for(int d=MINUS_ONE;++d < iDims;) cerr << iaMean[iNoOfDims*i+d] << ' ';
    cerr << endl;
  }
  for(int i=MINUS_ONE;++i < iNoOfClusters;) {
    for(int d=MINUS_ONE;++d < iDims;) cerr << daStdDev[iNoOfDims*i+d] << ' ';
    cerr << endl;
  }
#endif
}

void writeNoisyPoint(ofstream& fOut,int iNoOfDims,int iDims) {
  // for each dimension
  for(int d = MINUS_ONE;++d < iNoOfDims;) {
    int iRand = rand();
    if(d >= iDims) continue;
    fOut << iRand%RANGE;
    if(d != iDims - ONE) fOut << ' ';
  }
}

void writePoint(ofstream& fOut,int iNoOfDims,int iDims,int* iaMean,double *daStdDev,DistMode distMode) {
  double dRand;
  int iRand = 100;
  // for each dimension
  for(int d = MINUS_ONE;++d < iNoOfDims;) {
    dRand = genStdNormal();
    iRand=rand();
    if(d >= iDims) continue;
    if(iaMean[d]==MINUS_ONE) {
      fOut << iRand%RANGE;
    } else {
      switch(distMode) {
        case FULL_DIM:
          fOut << max(0,min(int(iaMean[d]+daStdDev[d]*dRand),RANGE));
          break;

        case HYPER_RECT_UNIFORM:
          fOut << max(0,min(RANGE,iaMean[d]-int(.5*daStdDev[d])+iRand%(int(daStdDev[d]))));
          break;

        case MVNORMAL:
          fOut << max(0,min(int(iaMean[d]+daStdDev[d]*dRand),RANGE));
          break;

        default:
          cerr << "ERROR: UNDEFINED DATA GENERATION MODE " << distMode << endl;
          exit(0);
      }
    }
    if(d != iDims - ONE) fOut << ' ';
  }
}

void writePt(ofstream& fOut,int iNoOfDims,int iDims,int* iaMean,double *daStdDev,int* iaCoord,DistMode distMode) {
  // choose a cluster
  int iRand = rand()%RANGE;
  int iClusterId = ZERO;
  while(iaCoord[iClusterId++] <= iRand);
  --iClusterId;
  writePoint(fOut,iNoOfDims,iDims,iaMean+iClusterId*iNoOfDims,daStdDev+iClusterId*iNoOfDims,distMode);
}

void parseArguments(int argc,char **argv) {

  // read args
  if(argc < NUM_ARGS) {
    cerr << "USAGE: dataGen\n";
    cerr << "\t-a<Max cluster size ratio   default="<<dMaxClusterSzRatio<<">\n";
    cerr << "\t-b<Min cluster size frac    default=" << dBias << ">\n";
    cerr << "\t-c<P(dim is bounded)        default=" << dProbConstrain << ">\n";
    cerr << "\t-d<Real Dims                default=" << iRealDims << ">\n";
    cerr << "\t-D<#(Dims)                  default=" << iNoOfDims << ">\n";
    cerr << "\t-e<generate subspaces only, no points>\n";
    cerr << "\t-m<mode                     default=" << iMode << ">\n";
    cerr << "\t-k<#(Clusters)              default=" << iNoOfClusters << ">\n";
    cerr << "\t-n<#(pts)                   default=" << iNoOfPts << ">\n";
    cerr << "\t-o<out file>\n";
    cerr << "\t-O<P(dim bounding overlaps) default=" << dProbOverlap << ">\n";
    cerr << "\t-r<noise                    default=" << dProbNoise << ">\n";
    cerr << "\t-R<range of variation       default=" << x << ">\n";
    cerr << "\t-s<#(datasets)              default=" << iNoOfDatasets << ">\n";
    cerr << "\t-w<width of clu in uni dist default=" << dWidth << ">\n";
    exit(0);
  }
  int c;
  while((c=getopt(argc,argv,"a:b:c:d:D:e:k:m:n:o:O:r:R:s:w:")) != -1) {
    switch(c) {
      case 'a':
        dMaxClusterSzRatio = atof(optarg);
        break;

      case 'b':
        dBias = atof(optarg);
        break;

      case 'c':
        dProbConstrain = atof(optarg);
        break;

      case 'd':
        iRealDims = atoi(optarg);
        break;

      case 'D':
        iNoOfDims = atoi(optarg);
        break;

      case 'e':
        bCentresOnly = true;
        break;

      case 'k':
        iNoOfClusters = atoi(optarg);
        break;

      case 'm':
	iMode = atoi(optarg);
        break;

      case 'n':
        iNoOfPts = atoi(optarg);
        break;

      case 'o':
	sprintf(caOutFile,"%s",optarg);
        break;

      case 'O':
	dProbOverlap = atof(optarg);
        break;

      case 'r':
        dProbNoise = atof(optarg);
        break;

      case 'R':
        x = atof(optarg);
        break;

      case 's':
        iNoOfDatasets = atoi(optarg);
        break;

      case 'w':
        dWidth = atof(optarg);
        break;

      default:
	cerr << "ERROR: OPTION -" << c << " NOT SUPPORTED\n";
	exit(0);
    }
  }
}

int main(int argc,char **argv) {
  parseArguments(argc,argv);
  DistMode distMode = DistMode(iMode);

  // test validity of parameters
  srand(20000);
  if(iRealDims > iNoOfDims) {
    cerr << " ERROR: INCORRECT DIMENSIONS : NoOfDims=" << iNoOfDims << " RealDims=" << iRealDims << endl;
    exit(0);
  }

  int k = (iNoOfDatasets > 1)? int((1.0+x)*iNoOfClusters):iNoOfClusters;
  int iaMean[k*iNoOfDims],iaMeanTmp[iNoOfDims];
  double daStdDev[k*iNoOfDims],daStdDevTmp[iNoOfDims];
  for(int i=ZERO;i<k*iNoOfDims;++i) {
    iaMean[i]=MINUS_ONE;
    daStdDev[i]=dStdDev;
  }
  createClusterParams(iaMean,daStdDev,iNoOfDims,k,distMode,iRealDims,dProbConstrain,dProbOverlap);
  for(int iDatasetId = ZERO;iDatasetId < iNoOfDatasets;++iDatasetId) {
    // open file to write to
    char caFile[1000],caOut[1000],caCentre[1000];
    memset(caOut,'\0',1000);
    if(iNoOfDatasets>ONE) {
      sprintf(caFile,"%.*ss%.2d%s",strlen(caOutFile)-3,caOutFile,iDatasetId,caOutFile+strlen(caOutFile)-3);
    } else sprintf(caFile,"%s",caOutFile);
    ofstream fHa(caFile,ios::out);
    if(fHa.bad()) {
      cerr << " ERROR: OPENING FILE " << caFile << endl;
      exit(0);
    }
    strncpy(caOut,caFile,strlen(caFile)-3);
    ofstream fOut(caOut,ios::out);
    if(fOut.bad()) {
      cerr << " ERROR: OPENING FILE " << caOut << endl;
      exit(0);
    }
    sprintf(caCentre,"%s_out",caOut);
    ofstream fCentre(caCentre,ios::out);
    if(fCentre.bad()) {
      cerr << " ERROR: OPENING FILE " << caCentre << endl;
      exit(0);
    }
    // perturb number of points and interesting subspaces in dataset
    int n = int(iNoOfPts*((iNoOfDatasets > 1)?x+ (rand()%int(1000*(1.0-x)))/500.0:1));
    int k = int(iNoOfClusters*((iNoOfDatasets > 1)?x+ (rand()%int(1000*(1.0-x)))/500.0:1));
    int iDims = iRealDims;//int(iRealDims*((iNoOfDatasets > 1)?x+ (rand()%int(1000*(1.0-x)))/500.0:1));
    int iaCoord[k];
    fCentre << iDims << endl;
    if(!bCentresOnly) {
      fOut << n << ' ' << iDims << ' ';
      fHa << iDims << endl;
      createClusterSizes(iaCoord,dBias,k,dMaxClusterSzRatio,dProbNoise,fOut);
    }
    int iPtsWritten = ZERO,iTmp = ZERO;
    for(int iCluId=ZERO;iCluId<k;++iCluId) { // for each cluster
      // determine number of points ascribed to that cluster
      int iNoOfCluPts = int(n*double(iaCoord[iCluId]-iTmp)/RANGE);
      iTmp = iaCoord[iCluId];
      for(int d=ZERO;d<iNoOfDims;++d) {
	// perturb the centre and deviation of the planted subspace
	if(iaMean[iNoOfDims*iCluId+d] < ZERO) {
	  iaMeanTmp[d]=MINUS_ONE;
          daStdDevTmp[d]=dStdDev;
	  continue;
	}
	double x = 0.02;
	iaMeanTmp[d]=max(0,min(RANGE,int(((iNoOfDatasets > 1)?iaMean[iCluId*iNoOfDims+d]-x*RANGE + 2*(rand()%int(RANGE*x)):iaMean[iCluId*iNoOfDims+d]))));
	daStdDevTmp[d]=daStdDev[iCluId*iNoOfDims+d]*((iNoOfDatasets > 1)?(1.0-x) + (rand()%int(1000*x))/500.0:1);
      }
      for(int d=MINUS_ONE;++d < iDims;) fCentre << iaMeanTmp[d] << ' ';
      fCentre << endl;
#ifdef DEBUG
//  for(int d=MINUS_ONE;++d < iDims;) fCentre << daStdDevTmp[d] << ' ';
//  fCentre << endl;
#endif
      if(!bCentresOnly) {
        for(int j=ZERO;j<iNoOfCluPts;++j) { // for each point ascribed to it
          if((iCluId+j) != ZERO) fHa << endl;
          writePoint(fHa,iNoOfDims,iDims,iaMeanTmp,daStdDevTmp,distMode);
        }
        iPtsWritten += iNoOfCluPts;
      }
    }
    if(!bCentresOnly) {
      while(++iPtsWritten <= n) {
        fHa << endl;
        writeNoisyPoint(fHa,iNoOfDims,iDims);
      }
      fHa << endl;
      fHa.flush();
    }
  }
  return 0;
}