##
## Bayesian and frequentist approaches to hypothesis testing for MISO
##
from numpy import *
import os
import scipy
from scipy import stats
from scipy.stats import gaussian_kde
from decimal import Decimal
import misopy
from misopy.samples_utils import *
from misopy.credible_intervals import *
import misopy.misc_utils as misc_utils

class NullPeakedDensity:
    """
    A density peaked on the null hypothesis
    """
    def __init__(self, dataset):
	self.dataset = dataset

    def evaluate(self, point):
	if point[0] == 0:
	    return inf
	else:
	    return 0

class gaussian_kde_set_covariance(stats.gaussian_kde):
    '''
    from Anne Archibald in mailinglist:
    http://www.nabble.com/Width-of-the-gaussian-in-stats.kde.gaussian_kde---td19558924.html#a19558924
    '''
    def __init__(self, dataset, covariance):
        self.covariance = covariance
        scipy.stats.gaussian_kde.__init__(self, dataset)

    def _compute_covariance(self):
        self.inv_cov = np.linalg.inv(self.covariance)
        self._norm_factor = sqrt(np.linalg.det(2*np.pi*self.covariance)) * self.n
        
class gaussian_kde_covfact(stats.gaussian_kde):
    def __init__(self, dataset, covfact = 'scotts'):
        self.covfact = covfact
        scipy.stats.gaussian_kde.__init__(self, dataset)

    def _compute_covariance_(self):
        '''not used'''
        self.inv_cov = np.linalg.inv(self.covariance)
        self._norm_factor = sqrt(np.linalg.det(2*np.pi*self.covariance)) * self.n

    def covariance_factor(self):
        if self.covfact in ['sc', 'scotts']:
            return self.scotts_factor()
        if self.covfact in ['si', 'silverman']:
            return self.silverman_factor()
        elif self.covfact:
            return float(self.covfact)
        else:
            raise ValueError, \
                'covariance factor has to be scotts, silv
erman or a number'

    def reset_covfact(self, covfact):
        self.covfact = covfact
        self.covariance_factor()
        self._compute_covariance()
	
def compute_prior_proportion_diff(num_samples):
    """
    Compute the posterior on the difference between two independent proportions (from two
    distinct conditions.)

    The model assumed here is:

    p_c1 ~ Beta(1, 1)
    p_c2 ~ Beta(1, 1)

    We compute P(delta), wher
e delta = p1_c1 - p2_c2, and return a set of num_samples-many samples.
    """
    samples = []
    for n in range(num_samples):
	# sample probabilities from uniform prior
	prob_c1 = random.beta(1, 1)
	prob_c2 = random.beta(1, 1)
	delta = prob_c1 - prob_c2
	samples.append(delta)
    return array(samples)


def compute_delta_densities(samples1_results,
                            samples2_results,
                            diff_range,
                            event_name="",
                            sample1_label="",
                            sample2_label="",
                            smoothing_param=0.3):
    """
    Compute the Gaussian kernel density fitted distributions
    over delta for the two sets of posterior samples filenames given.
    Returns the posterior density and prior density as well, assuming
    a uniform prior over the Psi of the samples
    in the two conditions.
    """
    densities = {}
    # Compute analytic prior density
    prior_density_fn = lambda x: 1 + x if x <= 0 else 1 - x
    analytic_prior_density = map(prior_density_fn, diff_range)
    posterior_samples1 = samples1_results[0]
    posterior_samples2 = samples2_results[0]

    num_samples, num_isoforms = shape(posterior_samples1)

    # Extract isoforms header information
    header = samples1_results[1]
    isoforms_field = get_isoforms_from_header(header[0])

    # Extract counts from headers
    sample1_counts_info = samples1_results[5]
    sample2_counts_info = samples2_results[5]

    # Record isoform information and counts
    densities['isoforms'] = isoforms_field
    densities['sample1_counts'] = sample1_counts_info
    densities['sample2_counts'] = sample2_counts_info

    # Set prior density function
    densities['prior_density'] = prior_density_fn

    # Posterior samples from MCMC
    densities['samples1'] = posterior_samples1
    densities['samples2'] = posterior_samples2
    
    # Collection of posterior densities (only 1 in two-isoform case)
    densities['posterior_density'] = []

    # Collection of Bayes factors (only 1 in two-isoform case)
    densities['bayes_factor'] = []
    
    # For each isoform, compute its Bayes factor and delta posterior
    warning_outputted = False
    
    for iso_num in range(num_isoforms):
        ##
        ## TODO: If distinct number of samples are used
        ## in the two samples compared, take the minimum
        ## number of samples
        ##
        posterior_diff = posterior_samples1[:, iso_num] - \
                         posterior_samples2[:, iso_num]

        # If the average difference is 0, don't try to fit a kernel to it
        mean_abs_posterior_diff = mean(abs(posterior_diff))

        # If the posterior differences are all identical, the sampler
        # was probably unable to explore the space
        all_same_diff = all(posterior_diff - posterior_diff[0] == 0)

        if all_same_diff and not warning_outputted:
            print "Warning: Event %s was not sampled properly in %s or %s" \
                  %(event_name,
                    sample1_label,
                    sample2_label)
            warning_outputted = True

        if mean_abs_posterior_diff <= .009 or all_same_diff:
            posterior_density = NullPeakedDensity(posterior_diff)
        else:
            # Smoothing by fitting Gaussian kernel density estimator
            posterior_density = gaussian_kde_covfact(posterior_diff,
                                                     smoothing_param)

        bayes_factor, diff_prior, diff_posterior = \
            compute_bayes_factor(analytic_prior_density,
                                 posterior_density)
        # Edge case: if improperly samples, maintain Bayes factors
        # as list to be consistent with multi-isoform case
        densities['bayes_factor'].append(bayes_factor)
        densities['posterior_density'].append(posterior_density)
    
    return densities


def output_samples_comparison(sample1_dir, sample2_dir, output_dir,
                              alpha=.95,
                              sample_labels=None,
                              use_compressed=None):
    """
    Compute the bayes factors, posterior means, and other statistics
    between the two samples and output them to a directory.

    Expects two directories with samples from a MISO run, where corresponding
    events in the two samples' directories begin with the same event name.
    """
    print "Given output dir: %s" %(output_dir)
    print "Retrieving MISO files in sample directories..."
    sample1_obj = MISOSamples(sample1_dir,
                              use_compressed=use_compressed)
    sample2_obj = MISOSamples(sample2_dir,
                              use_compressed=use_compressed)
    print "Computing sample comparison between %s and %s..." %(sample1_dir,
                                                               sample2_dir)
    print "  - No. of events in %s: %d" %(sample1_dir, sample1_obj.num_events)
    print "  - No. of events in %s: %d" %(sample2_dir, sample2_obj.num_events)
    # Output header for Bayes factor file
    if sample_labels is None:
        # Use directory names as sample labels
        sample1_label = os.path.basename(os.path.normpath(sample1_dir))
        sample2_label = os.path.basename(os.path.normpath(sample2_dir))
    else:
        # If we're given sample labels, use them
        sample1_label, sample2_label = sample_labels
        print "Using user-given sample labels (sample1 = %s, sample2 = %s)" \
              %(sample1_label, sample2_label)
    output_dir = os.path.join(output_dir, "%s_vs_%s" %(sample1_label,
                                                       sample2_label))
    print "Creating comparisons parent directory: %s" %(output_dir)
    # Create parent directory for comparison
    misc_utils.make_dir(output_dir)
	
    # Create d