CDFSampler.java

/*******************************************************************************
*   Copyright 2013 Analog Devices, Inc.
*
*   Licensed under the Apache License, Version 2.0 (the "License");
*   you may not use this file except in compliance with the License.
*   You may obtain a copy of the License at
*
*       http://www.apache.org/licenses/LICENSE-2.0
*
*   Unless required by applicable law or agreed to in writing, software
*   distributed under the License is distributed on an "AS IS" BASIS,
*   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*   See the License for the specific language governing permissions and
*   limitations under the License.
********************************************************************************/

package com.analog.lyric.dimple.solvers.gibbs.samplers.generic;

import static com.analog.lyric.dimple.environment.DimpleEnvironment.*;

import org.apache.commons.math3.random.RandomGenerator;

import com.analog.lyric.collect.ArrayUtil;
import com.analog.lyric.dimple.exceptions.DimpleException;
import com.analog.lyric.dimple.model.domains.DiscreteDomain;
import com.analog.lyric.dimple.model.domains.Domain;
import com.analog.lyric.dimple.model.values.DiscreteValue;
import com.analog.lyric.math.Utilities;

public class CDFSampler extends AbstractGenericSampler implements IDiscreteDirectSampler
{
	protected double[] _samplerScratch = ArrayUtil.EMPTY_DOUBLE_ARRAY;
	protected int _lengthRoundedUp = 0;
	protected int _length = 0;

	@Override
	public void initialize(Domain variableDomain)
	{
		int length = ((DiscreteDomain)variableDomain).size();
		_length = length;
	    _lengthRoundedUp = Utilities.nextPow2(length);
	    _samplerScratch = new double[_lengthRoundedUp];
	}
	
	@Override
	public void nextSample(DiscreteValue sampleValue, double[] energy, double minEnergy, IDiscreteSamplerClient samplerClient)
	{
		final RandomGenerator rand = activeRandom();
		final int length = sampleValue.getDomain().size(); //energy may be longer than domain size
		int sampleIndex;

		// Special-case lengths 2, 3, and 4 for speed
		switch (length)
		{
		case 2:
		{
			sampleIndex = (rand.nextDouble() * (1 + Math.exp(energy[1]-energy[0])) > 1) ? 0 : 1;
			break;
		}
		case 3:
		{
			final double cumulative1 = Math.exp(minEnergy-energy[0]);
			final double cumulative2 = cumulative1 + Math.exp(minEnergy-energy[1]);
			final double sum = cumulative2 + Math.exp(minEnergy-energy[2]);
			final double randomValue = sum * rand.nextDouble();
			sampleIndex = (randomValue > cumulative2) ? 2 : (randomValue > cumulative1) ? 1 : 0;
			break;
		}
		case 4:
		{
			final double cumulative1 = Math.exp(minEnergy-energy[0]);
			final double cumulative2 = cumulative1 + Math.exp(minEnergy-energy[1]);
			final double cumulative3 = cumulative2 + Math.exp(minEnergy-energy[2]);
			final double sum = cumulative3 + Math.exp(minEnergy-energy[3]);
			final double randomValue = sum * rand.nextDouble();
			sampleIndex = (randomValue > cumulative2) ? ((randomValue > cumulative3) ? 3 : 2) : ((randomValue > cumulative1) ? 1 : 0);
			break;
		}
		default:	// For all other lengths
		{
			// Calculate cumulative conditional probability (unnormalized)
			double sum = 0;
			final double[] samplerScratch = _samplerScratch;
			samplerScratch[0] = 0;
			for (int m = 1; m < length; m++)
			{
				sum += expApprox(minEnergy-energy[m-1]);
				samplerScratch[m] = sum;
			}
			sum += expApprox(minEnergy-energy[length-1]);
			for (int m = length; m < _lengthRoundedUp; m++)
				samplerScratch[m] = Double.POSITIVE_INFINITY;

			final int half = _lengthRoundedUp >> 1;
			while (true)
			{
				// Sample from the distribution using a binary search.
				final double randomValue = sum * rand.nextDouble();
				sampleIndex = 0;
				for (int bitValue = half; bitValue > 0; bitValue >>= 1)
				{
					final int testIndex = sampleIndex | bitValue;
					if (randomValue > samplerScratch[testIndex]) sampleIndex = testIndex;
				}

				// Rejection sampling, since the approximation of the exponential function is so coarse
				final double logp = minEnergy-energy[sampleIndex];
				if (Double.isNaN(logp))
					throw new DimpleException("The energy for all values of this variable is infinite. This may indicate a state inconsistent with the model.");
				if (rand.nextDouble()*expApprox(logp) <= Math.exp(logp)) break;
			}
		}
		}

		samplerClient.setNextSampleIndex(sampleIndex);
	}
	
	
	// This is an approximation to the exponential function; inputs must be non-positive
	// To facilitate subsequent rejection sampling, the error versus the correct exponential function needs to be always positive
	// This is true except for very large negative inputs, for values just as the output approaches zero
	// To ensure rejection is never in an infinite loop, this must reach 0 for large negative inputs before the Math.exp function does
	public final static double expApprox(double value)
	{
		// Convert input to base2 log, then convert integer part into IEEE754 exponent
		final long expValue = (long)((int)(1512775.395195186 * value) + 0x3FF00000) << 32;	// 1512775.395195186 = 2^20/log(2)
		return Double.longBitsToDouble(expValue & ~(expValue >> 63));	// Clip result if negative and convert to a double
	}
}