# Path: picrust/predict_traits.py
# def variance_of_weighted_mean(weights, variances, per_sample_axis=1):
#     """Calculate the standard deviation of a weighted average
# 
#     weights -- a numpy array of floats representing the weights
#     in the weighted average. Can be 1D or 2D (see per_sample_axis)
# 
#     variances -- a numpy array of variances for each observation (in
#     the same order as variance
# 
#     axis_for_samples -- which axis to sum down to stay within samples.
#     this value is ignored unless the array is 2d.
# 
#     This can be a little confusing, so consider this example:
#     If this is a 1D array, it is just handled asa single sample.
# 
#     If it is a 2D array, then multiple weighted averages will be calculated
#     for each column.   Lets say we set axis_for_samples equal to 1.
#     (axis = axis_for_samples).
# 
#     Now, if we wanted to simultaneously
#     calculate the variance of the weighted average for:
# 
#     case1: weights = [0.5,0.5], variances = [1.0,2.0]
#     case2: weights = [0.1,0.9], variances = [1000.0,2000.0]
# 
#     We would supply the following numpy arrays:
#     weights = [[0.5,0.5],[1.0,2.0]]
#     variances = [[1.0,2.0],[1000.0,2000.0]]
# 
#     Notes:
#     Formula from http://en.wikipedia.org/wiki/Weighted_mean
# 
#     The variance of the weighted mean is:
#     sqrt(sum(stdev_squared * weight_squared))
# 
#     Variance is standard deviation squared, so we use that instead
#     """
#     if variances.shape != weights.shape:
#         raise ValueError("Shape of variances (%s) doesn't match shape of weights (%s)" %(str(variances.shape),str(weights.shape)))
#     if len(variances.shape) > 2:
#         raise ValueError("Must supply a 1D or 2D array for variance.  Shape of supplied array was: %s" % str(variancs.shape))
#     if len(variances.shape) == 1:
#         per_sample_axis = None
#     #Ensure weights are normalized
# 
#     if per_sample_axis is not None:
#         normalized_weights = apply_along_axis(normalize_axis, per_sample_axis, weights)
#     else:
#         normalized_weights = normalize_axis(weights)
#     squared_weights = weights**2
# 
#     return  sqrt(sum(squared_weights*variances, axis=per_sample_axis))
# 
# def calc_confidence_interval_95(predictions,variances,round_CI=True,\
#         min_val=0.0,max_val=None):
#     """Calc the 95% confidence interval given predictions and variances"""
#     stdev = sqrt(variances)
#     pred = predictions
#     CI_95 =  1.96*stdev
#     lower_95_CI = pred - CI_95
#     upper_95_CI = pred + CI_95
#     if round_CI:
#         lower_95_CI = around(lower_95_CI)
#         upper_95_CI = around(upper_95_CI)
#     if min_val is not None:
#         lower_95_CI = numpy_max(min_val,lower_95_CI)
#         upper_95_CI = numpy_max(min_val,upper_95_CI)
# 
#     if max_val is not None:
#         lower_95_CI = numpy_min(max_val,lower_95_CI)
#         upper_95_CI = numpy_min(max_val,upper_95_CI)
# 
#     return lower_95_CI,upper_95_CI


from numpy import abs,compress, dot, array, around, asarray,empty,zeros, sum as numpy_sum,sqrt,apply_along_axis
from biom.table import Table
from biom.parse import parse_biom_table, get_axis_indices, direct_slice_data, direct_parse_key
from picrust.predict_traits import variance_of_weighted_mean,calc_confidence_interval_95
    metagenome_variance_data = None
    if verbose:
        print "Calculating the variance of the estimated metagenome for %i OTUs." %len(overlapping_otus)
    for otu_id in overlapping_otus:
        otu_across_samples = otu_table.data(otu_id, axis='observation')
        otu_across_genes = genome_table.data(otu_id)
        otu_variance_across_genes = gene_variances.data(otu_id)
        otu_contrib_to_metagenome = array([o*otu_across_genes
                                           for o in otu_across_samples])
        var_otu_contrib_to_metagenome = \
            array([scaled_variance(otu_variance_across_genes,o)
                   for o in otu_across_samples])

        if metagenome_data is None:
            metagenome_data = otu_contrib_to_metagenome
            metagenome_variance_data = var_otu_contrib_to_metagenome
        else:
            metagenome_data += otu_contrib_to_metagenome
            metagenome_variance_data = variance_of_sum(metagenome_variance_data,var_otu_contrib_to_metagenome)

    data_result = metagenome_data.T
    variance_result = metagenome_variance_data.T

    if whole_round:
        #Round counts to nearest whole numbers
        data_result = around(data_result)

    if verbose:
        print "Calculating metagenomic confidence intervals from variance."


# Path: picrust/predict_traits.py
# def variance_of_weighted_mean(weights, variances, per_sample_axis=1):
#     """Calculate the standard deviation of a weighted average
# 
#     weights -- a numpy array of floats representing the weights
#     in the weighted average. Can be 1D or 2D (see per_sample_axis)
# 
#     variances -- a numpy array of variances for each observation (in
#     the same order as variance
# 
#     axis_for_samples -- which axis to sum down to stay within samples.
#     this value is ignored unless the array is 2d.
# 
#     This can be a little confusing, so consider this example:
#     If this is a 1D array, it is just handled asa single sample.
# 
#     If it is a 2D array, then multiple weighted averages will be calculated
#     for each column.   Lets say we set axis_for_samples equal to 1.
#     (axis = axis_for_samples).
# 
#     Now, if we wanted to simultaneously
#     calculate the variance of the weighted average for:
# 
#     case1: weights = [0.5,0.5], variances = [1.0,2.0]
#     case2: weights = [0.1,0.9], variances = [1000.0,2000.0]
# 
#     We would supply the following numpy arrays:
#     weights = [[0.5,0.5],[1.0,2.0]]
#     variances = [[1.0,2.0],[1000.0,2000.0]]
# 
#     Notes:
#     Formula from http://en.wikipedia.org/wiki/Weighted_mean
# 
#     The variance of the weighted mean is:
#     sqrt(sum(stdev_squared * weight_squared))
# 
#     Variance is standard deviation squared, so we use that instead
#     """
#     if variances.shape != weights.shape:
#         raise ValueError("Shape of variances (%s) doesn't match shape of weights (%s)" %(str(variances.shape),str(weights.shape)))
#     if len(variances.shape) > 2:
#         raise ValueError("Must supply a 1D or 2D array for variance.  Shape of supplied array was: %s" % str(variancs.shape))
#     if len(variances.shape) == 1:
#         per_sample_axis = None
#     #Ensure weights are normalized
# 
#     if per_sample_axis is not None:
#         normalized_weights = apply_along_axis(normalize_axis, per_sample_axis, weights)
#     else:
#         normalized_weights = normalize_axis(weights)
#     squared_weights = weights**2
# 
#     return  sqrt(sum(squared_weights*variances, axis=per_sample_axis))
# 
# def calc_confidence_interval_95(predictions,variances,round_CI=True,\
#         min_val=0.0,max_val=None):
#     """Calc the 95% confidence interval given predictions and variances"""
#     stdev = sqrt(variances)
#     pred = predictions
#     CI_95 =  1.96*stdev
#     lower_95_CI = pred - CI_95
#     upper_95_CI = pred + CI_95
#     if round_CI:
#         lower_95_CI = around(lower_95_CI)
#         upper_95_CI = around(upper_95_CI)
#     if min_val is not None:
#         lower_95_CI = numpy_max(min_val,lower_95_CI)
#         upper_95_CI = numpy_max(min_val,upper_95_CI)
# 
#     if max_val is not None:
#         lower_95_CI = numpy_min(max_val,lower_95_CI)
#         upper_95_CI = numpy_min(max_val,upper_95_CI)
# 
#     return lower_95_CI,upper_95_CI


# Path: picrust/predict_metagenomes.py
from numpy import abs,compress, dot, array, around, asarray,empty,zeros, sum as numpy_sum,sqrt,apply_along_axis
from biom.table import Table
from biom.parse import parse_biom_table, get_axis_indices, direct_slice_data, direct_parse_key
from picrust.predict_traits import variance_of_weighted_mean,calc_confidence_interval_95

    metagenome_variance_data = None
    if verbose:
        print "Calculating the variance of the estimated metagenome for %i OTUs." %len(overlapping_otus)
    for otu_id in overlapping_otus:
        otu_across_samples = otu_table.data(otu_id, axis='observation')
        otu_across_genes = genome_table.data(otu_id)
        otu_variance_across_genes = gene_variances.data(otu_id)
        otu_contrib_to_metagenome = array([o*otu_across_genes
                                           for o in otu_across_samples])
        var_otu_contrib_to_metagenome = \
            array([scaled_variance(otu_variance_across_genes,o)
                   for o in otu_across_samples])

        if metagenome_data is None:
            metagenome_data = otu_contrib_to_metagenome
            metagenome_variance_data = var_otu_contrib_to_metagenome
        else:
            metagenome_data += otu_contrib_to_metagenome
            metagenome_variance_data = variance_of_sum(metagenome_variance_data,var_otu_contrib_to_metagenome)

    data_result = metagenome_data.T
    variance_result = metagenome_variance_data.T

    if whole_round:
        #Round counts to nearest whole numbers
        data_result = around(data_result)

    if verbose:
        print "Calculating metagenomic confidence intervals from variance."


    lower_95_CI,upper_95_CI=calc_confidence_interval_95(data_result,variance_result,\