//==============================================================================
//	
//	Copyright (c) 2002-
//	Authors:
//	* Dave Parker <david.parker@comlab.ox.ac.uk> (University of Oxford)
//	* Vincent Nimal <vincent.nimal@comlab.ox.ac.uk> (University of Oxford)
//	* Marcin Copik <mcopik@gmail.com>
//	
//------------------------------------------------------------------------------
//	
//	This file is part of PRISM.
//	
//	PRISM 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 2 of the License, or
//	(at your option) any later version.
//	
//	PRISM 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 PRISM; if not, write to the Free Software Foundation,
//	Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//	
//==============================================================================

package simulator.sampler;

import simulator.*;
import prism.PrismException;
import prism.PrismLangException;

/**
 * Samplers for properties that associate a simulation path with a real (double) value.
 */
public abstract class SamplerDouble extends Sampler
{
	/** Value of current path */
	protected double value;
	
	// Stats over all paths
	
	/** Number of samples so far */
	protected int numSamples;
	/** Sum of values (used to compute mean) */
	protected double valueSum;
	/** Correction value used when tracking variance */
	protected double correctionTerm;
	/** Sum of values, each shifted by the correction term (see below) */
	protected double valueSumShifted;
	/** Sum of squares of values, each shifted by the correction term (see below) */
	protected double valueSumShiftedSq;

	/**
	 * NB: In order to improve the numerical stability for the computation of variance,
	 * we stored a shifted version of the mean, using a correction term that is an
	 * estimate of the mean (for which, we just use the first value sampled). 
	 * Source: "Algorithms for Computing the Sample Variance: Analysis and Recommendations", T.F. Chan et al.
	 */
	
	@Override
	public void reset()
	{
		valueKnown = false;
		value = 0.0;
	}

	@Override
	public void resetStats()
	{
		valueSum = 0.0;
		valueSumShifted = 0.0;
		valueSumShiftedSq = 0.0;
		numSamples = 0;
	}

	@Override
	public abstract boolean update(Path path, TransitionList transList) throws PrismLangException;

	@Override
	public void updateStats()
	{
		if (numSamples == 0)
			correctionTerm = value;
		valueSum += value;
		valueSumShifted += value - correctionTerm;
		valueSumShiftedSq += Math.pow(value - correctionTerm, 2);
		numSamples++;
	}

	@Override
	public Object getCurrentValue()
	{
		return new Double(value);
	}

	@Override
	public double getMeanValue()
	{
		return valueSum / numSamples;
	}

	@Override
	public double getVariance()
	{
		// Return estimator to the variance
		if (numSamples <= 1) {
			return 0.0;
		} else {
			double meanShifted = valueSumShifted / numSamples;
			return (valueSumShiftedSq - numSamples * meanShifted * meanShifted) / (numSamples - 1.0);
		}
		
		// Note:
		// * As explained at the top of this class, variance is computed using
		//   a shifted version of the mean to improve numerical stability
		// * We divide by numSamples -1. An alternative would be to use the empirical mean, i.e. dividing by numSamples. 
		//   This is not equivalent (or unbiased) but, asymptotically, is the same.
		//   See p.24 of: Vincent Nimal, "Statistical Approaches for Probabilistic Model Checking", MSc Thesis
	}

	@Override
	public double getLikelihoodRatio(double p1, double p0) throws PrismException
	{
		// For details, see Sec 6.3 of: Vincent Nimal, "Statistical Approaches for Probabilistic Model Checking", MSc Thesis
		// (in which mu1=p1 and mu0=p0)
		if (numSamples <= 1)
			return 0.0;
		if (valueSumShiftedSq == 0)
			throw new PrismException("Cannot compute likelihood ratio with null variance");
		// Compute maximum likelihood estimator of variance
		double MLE = valueSumShiftedSq / numSamples - (valueSumShifted * valueSumShifted) / numSamples / numSamples;
		double lr = (-1 / (2 * MLE)) * (numSamples * (p1 * p1 - p0 * p0) - 2 * valueSum * (p1 - p0));
		if (Double.isNaN(lr)) {
			throw new PrismException("Error computing likelihood ratio");
		}
		return Math.exp(lr);
	}
}
