/**
*   Make up the Q and A arrays for Walker's Alias method for sampling
* from a discrete distribution with M terms.
*   Methodology: sort the distribution from low to high, move
* probability from the highest to the lowest to bring the latter up
* to 1/M; remove that cell and repeat with M-1 cells. 
*/
import java.util.*;
import static java.lang.System.err;

public class
walkerAlias
{
  static Random rndm;
  static long seed;
  static int nTrials = 20000;

//ggg
  double[] Q = null;
  int[] A = null;
  
//ccc===
  public
  walkerAlias(double[] targ)
  {
    int M = targ.length;
    Q = new double[M];
    A = new int[M];
    ArrayList<double[]> rTarg = new ArrayList(); //remaining to allocate
    for( int i=0; i<M; ++i)
      rTarg.add(new double[] {targ[i],(double)i});
    int[] rank = null;

    int j,k, p, q;
    double td; int ti;
    double leftover = 0;
    double[] dnrj = null, dnrk = null; //dataNindex
    for( int i=0; i<M; ++i)
    {
        //sort rTarg smallest to largest
      rank = new int[rTarg.size()];
      sort(rTarg,rank);
   
      j = rank[0]; k = rank[-1+rTarg.size()];
      dnrj = rTarg.get(j); dnrk = rTarg.get(k);
//err.format("least: %.3f(%d), most: %.3f(%d)\n",
//dnrj[0],(int)dnrj[1],dnrk[0],(int)dnrk[1]);
    
      p = (int)dnrj[1]; q = (int)dnrk[1];
      Q[p] = dnrj[0]*M;
      leftover = (1-Q[p])/M;
      A[p] = q;
      dnrk[0] -= leftover; rTarg.set(k,dnrk);
//err.format("Q[%d]= %.3f, A[%d]= %d, lo= %.3f\n",
//(int)dnrj[1],Q[p],p,A[p],leftover);

        //take index j out of rTarg, 
      rTarg.remove(j);
//for( int o=0; o<rTarg.size(); ++o ){
//dnrj = rTarg.get(o);
//err.format("  %.3f(%d) ",dnrj[0],p);}
//err.println("");
    }

for( int i=0; i<M; ++i)
err.format("  Q[%d]= %.3f",i,Q[i]);
err.println("");
for( int i=0; i<M; ++i)
err.format("  A[%d]= %d",i,A[i]);
err.println("");

  }
//mmm===
  public static void
  main(String[] args)
  {
    double[] targ = new double[] {
//0.1, 0.25, 0.15, 0.35, 0.15};
0.02, 0.04, 0.08, 0.15, 0.2, 0.21, 0.16, 0.09, 0.05}; //sum = 1.0

    walkerAlias wa = new walkerAlias(targ);

      //test it
    seed = System.currentTimeMillis(); //19;
err.println("seed= "+seed);
    rndm = new Random(seed);

    int X, Y;
    double U, Z;
    int M = targ.length;
    int[] histo = new int[M];

    for( int i=0; i<nTrials; ++i)
    {
      U = rndm.nextDouble()*M;
      Y = (int)U; Z = U-Y;
      if( Z < wa.Q[Y]) X=Y; else X=wa.A[Y];
      ++histo[X];
//err.println("Y= "+Y+", Z= "+Z+", X= "+X);
    }

    System.out.println("Empirical");
    for( int i=0; i<M; ++i)
      System.out.format("  %.2f",histo[i]/(double)nTrials);
    System.out.println("");
    System.out.println("Target");
    for( int i=0; i<M; ++i)
      System.out.format("  %.2f",targ[i]);
    System.out.println("");

  }
//===
   /**
   * first coor = data, 2nd = index
   */
  static void
  sort(ArrayList<double[]> data, int rank[])
  {
    int gap, len=data.size(); boolean done;
    int temp=0;
    for( int i=0; i<len; i++) rank[i] = i;

    double[] dnr1 = null, dnr2 = null;
    gap=len/2;
    while( gap>=1 )
    {
      do {
        done=true;
        for( int i=0; i<len-gap; i++)
        {
          dnr1 = data.get(rank[i]);
          dnr2 = data.get(rank[i+gap]);
          if( dnr1[0] > dnr2[0] )
          {
            temp = rank[i];
            rank[i] = rank[i+gap];
            rank[i+gap] = temp;
            done = false;
          } //end of if
        } //end of for
      } while(!done); //end of do.while
      gap = gap/2;
    }  //end of while
//err.println("Sorted data:");
//for( int k=0; k<data.length; k++) {
//err.println(data.get(rank[k])+" "); }
  }
}