/*
* #%L
* org.gitools.analysis
* %%
* Copyright (C) 2013 - 2014 Universitat Pompeu Fabra - Biomedical Genomics group
* %%
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program. If not, see
* <http://www.gnu.org/licenses/gpl-3.0.html>.
* #L%
*/
package org.gitools.analysis.stats.test;
/*
* Copyright (c) 2003-2005 The BISON Project
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License version 2 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*/
import java.util.NoSuchElementException;
import java.util.Random;
/**
* This class provides a weighted random permutation of indexes.
* Useful for weighted random sampling without replacement.
* The next sample is taken according to the weights given as a parameter
* to {@link #reset(int)}.
* The weights work as follows.
* The first sample is drawn according to the probability distribution
* defined by the (normalized) weights.
* After this the remaining k-1 elements and the associated k-1
* (re-normalized) weights
* define a new probability distribution, according to which the 2nd element
* is drawn, and so on.
*/
public class WeightedRandPerm {
// ======================= private fields ============================
// ===================================================================
/** Holds the weights that are used to initialize the permutation */
private final double[] w;
/** Holds the sum of the weights until the given index, inclusive. */
private final double[] wsum;
private int[] buffer = null;
/** Working array for calculating the permutation */
private double[] weights = null;
private int len = 0;
private int pointer = 0;
private double sum = 0.0;
private final Random r;
// ======================= initialization ============================
// ===================================================================
/** Set the source of randomness to use and the weights. You need to call
* {@link #reset} to fully initialize the object.
* @param r source of randomness
* @param weights The array that holds the weights for the calculation of the
* permutation. The length of the array will be an upper bound on the
* parameter {@link #reset} accepts. If {@link #reset} is called with a
* parameter less than the length of weights, the prefix of the same length
* is used.
* The vector elements must be positive, that is, zero is not accepted either.
*/
public WeightedRandPerm( Random r, double[] weights ) {
this.r=r;
w = weights.clone();
wsum = weights.clone();;
this.weights = new double[w.length];
buffer = new int[w.length];
for(int i=0; i<w.length; ++i)
{
if( w[i] <= 0.0 ) throw new IllegalArgumentException(
"weights should be positive: w["+i+"]="+w[i]);
}
for(int i=1; i<w.length; ++i) wsum[i]+=wsum[i-1];
}
// ======================= public methods ============================
// ===================================================================
/**
* It initiates a random weighted permutation of the integeres from 0 to k-1.
* It does not actually calculate the permutation.
* The permutation can be read using method {@link #next}.
* If the previous permutation was of the same length, it is more efficient.
* The weights set at construction time work as follows.
* The first sample is drawn according to the probability distribution
* defined by the (normalized) weights.
* After this the remaining k-1 elements and the associated k-1
* (re-normalized) weights
* define a new probability distribution, according to which the 2nd element
* is drawn, and so on.
* @param k the set is defined as 0,...,k-1
*/
public void reset(int k) {
if( k<0 || k>w.length )
throw new IllegalArgumentException(
"k should be non-negative and <= "+w.length);
pointer = k;
sum = wsum[k-1];
if( k != len )
{
// we need to initialize weights and buffer
for(int i=0; i<k; ++i)
{
weights[i]=w[i];
buffer[i]=i;
}
len=k;
}
}
// -------------------------------------------------------------------
/**
* The first sample is drawn according to the probability distribution
* defined by the (normalized) weights.
* After this the remaining k-1 elements and the associated k-1
* (re-normalized) weights
* define a new probability distribution, according to which the 2nd element
* is drawn, and so on.
* @see #reset
*/
public int next() {
if( pointer < 1 ) throw new NoSuchElementException();
double d = sum*r.nextDouble();
int i = pointer;
double tmp = weights[i-1];
while( tmp < d && i>1 ) tmp += weights[--i-1];
// now i-1 is the selected element, we shift it to next position
int a = buffer[i-1];
double b = weights[i-1];
buffer[i-1] = buffer[pointer-1];
weights[i-1] = weights[pointer-1];
buffer[pointer-1] = a;
weights[pointer-1] = b;
sum -= b;
return buffer[--pointer];
}
// -------------------------------------------------------------------
public boolean hasNext() { return pointer > 0; }
// -------------------------------------------------------------------
/*
public static void main( String pars[] ) throws Exception {
int k = pars.length;
double w[] = new double[k];
for(int i=0; i<k; ++i) w[i] = Double.parseDouble(pars[i]);
WeightedRandPerm rp = new WeightedRandPerm(new Random(),w);
rp.reset(k);
for(int i=0; i<1000; ++i)
{
if(i%2==0) rp.reset(k);
if(i%2==1) rp.reset(k-1);
while(rp.hasNext()) System.out.print(rp.next()+" ");
System.out.println();
}
System.out.println();
}
*/
}