<commit_before>from __future__ import division
import numpy as np
from scipy.stats import norm
from ..util import read_param_file
import common_args

# Perform Sobol Analysis on file of model results
# Returns a dictionary with keys 'S1', 'S1_conf', 'ST', and 'ST_conf'
# Where each entry is a list of size D (the number of parameters)
# Containing the indices in the same order as the parameter file
def analyze(pfile, output_file, column = 0, calc_second_order = True, num_resamples = 1000,
            delim = ' ', conf_level = 0.95, print_to_console=False):
    
    param_file = read_param_file(pfile)
    Y = np.loadtxt(output_file, delimiter=delim, usecols=(column,))
    D = param_file['num_vars']

    if calc_second_order and Y.size % (2*D + 2) == 0:
        N = int(Y.size / (2*D + 2))
    elif not calc_second_order and Y.size % (D + 2) == 0:
        N = int(Y.size / (D + 2))
    else: raise RuntimeError("""
        Incorrect number of samples in model output file. 
        Confirm that calc_second_order matches option used during sampling.""")  

    if conf_level < 0 or conf_level > 1: raise RuntimeError("Confidence level must be between 0-1.")
    
    A = np.empty(N)
    B = np.empty(N)
    AB = np.empty((N,D))
    BA = np.empty((N,D)) if calc_second_order else None
    step = 2*D+2 if calc_second_order else D+2

    A = Y[0:Y.size:step]
    B =  Y[(step-1):Y.size:step]
    for j in range(D):
        AB[:,j] = Y[(j+1):Y.size:step]
        if calc_second_order: BA[:,j] = Y[(j+1+D):Y.size:step]
    
    # First order (+conf.) and Total order (+conf.)
    keys = ('S1','S1_conf','ST','ST_conf')
    S = dict((k, np.empty(D)) for k in keys)
    if print_to_console: print "Parameter %s %s %s %s" % keys

    for j in range(D):
        S['S1'][j] = compute_first_order(A, AB[:,j], B, N)
        S['S1_conf'][j] = compute_first_order_confidence(A, AB[:,j], B, N, num_resamples, conf_level)
        S['ST'][j] = compute_total_order(A, AB[:,j], B, N)
        S['ST_conf'][j] = compute_total_order_confidence(A, AB[:,j], B, N, num_resamples, conf_level)
        
        if print_to_console:
            print "%s %f %f %f %f" % (param_file['names'][j], S['S1'][j], S['S1_conf'][j], S['ST'][j], S['ST_conf'][j])
    
    # Second order (+conf.)
    if calc_second_order:
        S['S2'] = np.empty((D,D))
        S['S2_conf'] = np.empty((D,D))
        if print_to_console: print "\nParameter_1 Parameter_2 S2 S2_conf"
        
        for j in range(D):
            for k in range(j+1, D):
                    
                S['S2'] = compute_second_order(A, BA[:,j], AB[:,k], AB[:,j], B, N)
                S['S2_conf'] = compute_second_order_confidence(A, BA[:,j], AB[:,k], AB[:,j], B, N, num_resamples, conf_level)
                
                if print_to_console:
                    print "%s %s %f %f" % (param_file['names'][j], param_file['names'][k], S['S2'], S['S2_conf'])                        
    
    return S            
        
def compute_first_order(a0, a1, a2, N):

    c = np.mean(a0)
    EY2 = np.mean((a0-c)*(a2-c))
    V = np.sum((a2-c)**2)/(N-1) - np.mean(a2-c)**2
    U = np.sum((a1-c)*(a2-c))/(N-1)

    return (U - EY2) / V

def compute_first_order_confidence(a0, a1, a2, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])
    
    for i in range(num_resamples):        
        r = np.random.randint(N, size=N)        
        s[i] = compute_first_order(a0[r], a1[r], a2[r], N)
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

def compute_total_order(a0, a1, a2, N):
    
    c = np.average(a0)
    EY2 = np.mean(a0-c)**2
    V = np.sum((a0-c)**2)/(N-1) - EY2
    U = np.sum((a0-c)*(a1-c))/(N-1)

    return (1 - (U-EY2) / V)

def compute_total_order_confidence(a0, a1, a2, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])

    if conf_level < 0 or conf_level > 1:    
        print "Error: Confidence level must be between 0-1."
        exit() 
    
    for i in range(num_resamples):
        r = np.random.randint(N, size=N)  
        s[i] = compute_total_order(a0[r], a1[r], a2[r], N)
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

def compute_second_order(a0, a1, a2, a3, a4, N):
    
    c = np.average(a0)    
    EY = np.mean((a0-c)*(a4-c))
    EY2 = np.mean((a1-c)*(a3-c))
    
    V = np.sum((a1-c)**2)/(N-1) - np.mean(a1-c)**2
    Vij = np.sum((a1-c)*(a2-c))/(N-1) - EY2
    Vi = np.sum((a2-c)*(a4-c))/(N-1) - EY
    Vj = np.sum((a3-c)*(a4-c))/(N-1) - EY2
    
    return (Vij - Vi - Vj) / V

def compute_second_order_confidence(a0, a1, a2, a3, a4, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])

    if conf_level < 0 or conf_level > 1:    
        print "Error: Confidence level must be between 0-1."
        exit() 
    
    for i in range(num_resamples):
        r = np.random.randint(N, size=N)
        s[i] = compute_second_order(a0[r], a1[r], a2[r], a3[r], a4[r], N)
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

if __name__ == "__main__":

    parser = common_args.create()
    parser.add_argument('--max-order', type=int, required=False, default=2, choices=[1, 2], help='Maximum order of sensitivity indices to calculate')
    parser.add_argument('-r', '--resamples', type=int, required=False, default=1000, help='Number of bootstrap resamples for Sobol confidence intervals')
    args = parser.parse_args()

    analyze(args.paramfile, args.model_output_file, args.column, (args.max_order == 2), num_resamples = args.resamples, delim = args.delimiter, print_to_console=True)
<commit_msg>first order estimator saltelli et al 2010<commit_after>from __future__ import division
import numpy as np
from scipy.stats import norm
from ..util import read_param_file
import common_args

# Perform Sobol Analysis on file of model results
# Returns a dictionary with keys 'S1', 'S1_conf', 'ST', and 'ST_conf'
# Where each entry is a list of size D (the number of parameters)
# Containing the indices in the same order as the parameter file
def analyze(pfile, output_file, column = 0, calc_second_order = True, num_resamples = 1000,
            delim = ' ', conf_level = 0.95, print_to_console=False):
    
    param_file = read_param_file(pfile)
    Y = np.loadtxt(output_file, delimiter=delim, usecols=(column,))
    D = param_file['num_vars']

    if calc_second_order and Y.size % (2*D + 2) == 0:
        N = int(Y.size / (2*D + 2))
    elif not calc_second_order and Y.size % (D + 2) == 0:
        N = int(Y.size / (D + 2))
    else: raise RuntimeError("""
        Incorrect number of samples in model output file. 
        Confirm that calc_second_order matches option used during sampling.""")  

    if conf_level < 0 or conf_level > 1: raise RuntimeError("Confidence level must be between 0-1.")
    
    A = np.empty(N)
    B = np.empty(N)
    AB = np.empty((N,D))
    BA = np.empty((N,D)) if calc_second_order else None
    step = 2*D+2 if calc_second_order else D+2

    A = Y[0:Y.size:step]
    B =  Y[(step-1):Y.size:step]
    for j in range(D):
        AB[:,j] = Y[(j+1):Y.size:step]
        if calc_second_order: BA[:,j] = Y[(j+1+D):Y.size:step]
    
    # First order (+conf.) and Total order (+conf.)
    keys = ('S1','S1_conf','ST','ST_conf')
    S = dict((k, np.empty(D)) for k in keys)
    if print_to_console: print "Parameter %s %s %s %s" % keys

    for j in range(D):
        S['S1'][j] = compute_first_order(A, AB[:,j], B)
        S['S1_conf'][j] = compute_first_order_confidence(A, AB[:,j], B, N, num_resamples, conf_level)
        S['ST'][j] = compute_total_order(A, AB[:,j], B, N)
        S['ST_conf'][j] = compute_total_order_confidence(A, AB[:,j], B, N, num_resamples, conf_level)
        
        if print_to_console:
            print "%s %f %f %f %f" % (param_file['names'][j], S['S1'][j], S['S1_conf'][j], S['ST'][j], S['ST_conf'][j])
    
    # Second order (+conf.)
    if calc_second_order:
        S['S2'] = np.empty((D,D))
        S['S2_conf'] = np.empty((D,D))
        if print_to_console: print "\nParameter_1 Parameter_2 S2 S2_conf"
        
        for j in range(D):
            for k in range(j+1, D):     
                S['S2'] = compute_second_order(A, BA[:,j], AB[:,k], AB[:,j], B, N)
                S['S2_conf'] = compute_second_order_confidence(A, BA[:,j], AB[:,k], AB[:,j], B, N, num_resamples, conf_level)
                
                if print_to_console:
                    print "%s %s %f %f" % (param_file['names'][j], param_file['names'][k], S['S2'], S['S2_conf'])                        
    
    return S            
        
def compute_first_order(A, AB, B):
    # First order estimator following Saltelli et al. 2010 CPC, normalized by sample variance
    return np.mean(B*(AB-A))/np.var(np.r_[A,B])

def compute_first_order_confidence(a0, a1, a2, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])
    
    for i in range(num_resamples):        
        r = np.random.randint(N, size=N)        
        s[i] = compute_first_order(a0[r], a1[r], a2[r])
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

def compute_total_order(a0, a1, a2, N):
    
    c = np.average(a0)
    EY2 = np.mean(a0-c)**2
    V = np.sum((a0-c)**2)/(N-1) - EY2
    U = np.sum((a0-c)*(a1-c))/(N-1)

    return (1 - (U-EY2) / V)

def compute_total_order_confidence(a0, a1, a2, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])

    if conf_level < 0 or conf_level > 1:    
        print "Error: Confidence level must be between 0-1."
        exit() 
    
    for i in range(num_resamples):
        r = np.random.randint(N, size=N)  
        s[i] = compute_total_order(a0[r], a1[r], a2[r], N)
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

def compute_second_order(a0, a1, a2, a3, a4, N):
    
    c = np.average(a0)    
    EY = np.mean((a0-c)*(a4-c))
    EY2 = np.mean((a1-c)*(a3-c))
    
    V = np.sum((a1-c)**2)/(N-1) - np.mean(a1-c)**2
    Vij = np.sum((a1-c)*(a2-c))/(N-1) - EY2
    Vi = np.sum((a2-c)*(a4-c))/(N-1) - EY
    Vj = np.sum((a3-c)*(a4-c))/(N-1) - EY2
    
    return (Vij - Vi - Vj) / V

def compute_second_order_confidence(a0, a1, a2, a3, a4, N, num_resamples, conf_level):
    
    s  = np.empty([num_resamples])

    if conf_level < 0 or conf_level > 1:    
        print "Error: Confidence level must be between 0-1."
        exit() 
    
    for i in range(num_resamples):
        r = np.random.randint(N, size=N)
        s[i] = compute_second_order(a0[r], a1[r], a2[r], a3[r], a4[r], N)
    
    return norm.ppf(0.5 + conf_level/2) * s.std(ddof=1)

if __name__ == "__main__":

    parser = common_args.create()
    parser.add_argument('--max-order', type=int, required=False, default=2, choices=[1, 2], help='Maximum order of sensitivity indices to calculate')
    parser.add_argument('-r', '--resamples', type=int, required=False, default=1000, help='Number of bootstrap resamples for Sobol confidence intervals')
    args = parser.parse_args()

    analyze(args.paramfile, args.model_output_file, args.column, (args.max_order == 2), num_resamples = args.resamples, delim = args.delimiter, print_to_console=True)

{'fraction_non_alphanumeric': 0.16117440841367223, 'fraction_numerical': 0.02646801051709027, 'mean_word_length': 2.4474320241691845, 'pattern_counts': {'":': 2, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6574572', 'n_tokens_mistral': 4573, 'n_tokens_neox': 4322, 'n_words': 1017}