//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// WAM = Walker's Alias Method for sampling discrete distribution
// NON-STATIC, OBJECT-ORIENTED VERSION.
//
// To compile: javac wamdemo.java      or "m wamdemo"
// To run:     java wamdemo
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
import java.io.*;
import java.util.Random;


//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// YLL 2005/05/07
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
class WAMsampler {  // N=number of sample values, ff=bar fractions, aa=aliases
  Random random; int N; double ff[]; int aa[];
  public WAMsampler (final int N, final double pp[]) { //pp=probability vector
    int i,j,jmin,jmax;
    double bmin,bmax,s,oon,tol = 1E-15d;
    double bb[] = new double[N];
    random = new Random();
    this.N = N;
    ff = new double[N];
    aa = new int[N];
    oon = 1.0d/N;
    //----- Verify that user-supplied pp[]'s sum to unity!
    s=0;
    for (i=0; i<N; i++) {s+=pp[i];}
    if (Math.abs(s-1) > tol) {System.out.println("not norm!");Runtime.getRuntime().exit(1);}
    //----- Set up arrays
    for (i=0; i<N; i++) {bb[i]=pp[i]-oon; aa[i]=i; ff[i]=1.0;}
    for (i=0; i<N; i++) {
      bmin=+1E10; jmin=-1; bmax=-1E10; jmax=-1;
      for (j=0; j<N; j++) {
	if (bmin>bb[j]) {bmin=bb[j]; jmin=j;}
	if (bmax<bb[j]) {bmax=bb[j]; jmax=j;}
      }
      if (bmax-bmin<tol) break;
      aa[jmin]=jmax; ff[jmin]=1+bmin*N; bb[jmin]=0; bb[jmax]=bmax+bmin;
    }
  }

  public int sample () {
    int n=random.nextInt(N); if (random.nextDouble()>ff[n]) n=aa[n]; return n;
  }
}


//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// TEST DRIVE
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
public class wamdemo { 
  public static void main (String[] args) {
    WAMsampler wamsampler;
    int N;
    double pp[]={.37, .04, .11, .06, .08, .19, .03, .01, .02, .09};
    N = 10;
    int i;
    wamsampler = new WAMsampler(N, pp);

    //------- Dump internals
    System.out.println ("\ni \tpp[i] \tff[i] \t1-f\taa[i]");
    for (i=0; i<N; i++) {
      System.out.println (i+"\t"+pp[i]+"\t"+
			  Math.round(100*wamsampler.ff[i])
			  +"%\t"+
			  Math.round(100*(1-wamsampler.ff[i]))
			  +"%\t"+wamsampler.aa[i]);
    }
    //---- Test drive and make histogram
    for (i=0; i<N; i++) pp[i]=0;      
    for (i=0; i<10000000; i++) {
      int n=wamsampler.sample();  //System.out.print (n + "");
      pp[n]++;
    }
    System.out.println ();
    for (i=0; i<N; i++) 
      System.out.println (i + "\t" + Math.round(pp[i]));
  }
}