<reponame>griffij/OtagoCosmo_publish
"""Sample age uncertainty distributions in order to get combined estimate
of age of surface"""

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
import matplotlib

#Rock Creek

data = np.genfromtxt('./data/test_ages.csv', delimiter=',', skip_header=1) [:,2:]
surface_labels = np.genfromtxt('./data/test_ages.csv', delimiter=',', skip_header=1, usecols=0, dtype=str)

# Get samples we want from each surface
# Edit the code below to add additional surfaces.
# Name must match surface labels in first column of data file
age_dict = {'S1':[]}#, 'T2':[], 'T3':[]}
sigma_dict = {'S1':[]}#, 'S2':[], 'S3':[]}
sample_dict = {'S1':[]}#, 'S2':[], 'S3':[]}
for i, label in enumerate(surface_labels):
    if label == 'S1':
        age_dict['S1'].append(data[i][0])
        sigma_dict['S1'].append(data[i][1])
#    elif label == 'S2':
#        age_dict['S2'].append(data[i][0])
#        sigma_dict['S2'].append(data[i][1]) 
#    elif label == 'S3':
#        age_dict['S3'].append(data[i][0])
#        sigma_dict['S3'].append(data[i][1])
    else:
        print('Sample from unknown surface')

# Check we've read it in correctly
print(age_dict)
print(sigma_dict)
for key, value in age_dict.items():
    for i, mu in enumerate(value):
        dnorm = norm(loc=mu, scale=sigma_dict[key][i])
        samples = dnorm.rvs(size=10000)
        sample_dict[key].append(samples)
    # Make an example plot
    if key=='S1':
        # Estimate bounds from data, may need to adjust manually (see line below)
        xvals = np.arange(min(value)-6*max(sigma_dict[key]), max(value)+6*max(sigma_dict[key]), 0.1)
#        xvals = np.arange(0, 40, 0.1)
        for i, mu in enumerate(value):
            dnorm = norm(loc=mu, scale=sigma_dict[key][i])
            plt.plot(xvals, dnorm.pdf(xvals), c='0.6')
            p1, = plt.fill(xvals, dnorm.pdf(xvals), facecolor='0.8', label='Individual sample') 
        mu_samples = np.mean(sample_dict[key])
        sig_samples = np.std(sample_dict[key])
        dnorm_samp = norm(loc=mu_samples, scale=sig_samples)
        plt.plot(xvals, dnorm_samp.pdf(xvals), c='0.3', zorder=11)
        p2, = plt.fill(xvals, dnorm_samp.pdf(xvals), facecolor='0.4',
                       zorder=10, label = 'Combined age')#, alpha=0.5)
        plt.xlabel('Age (ka)')
        plt.ylabel('Density')
        plt.legend(handles=[p1,p2])
        plt.savefig('plots/sampling_method_test_example.png', dpi=300)
        
print('\nTest data')
for key, value in sample_dict.items():
    mean = np.mean(value)
    sigma = np.std(value)
    print(key)
    print('Mean age (ka)', mean)
    print('Sigma (kyr)', sigma)

"""Sample age uncertainty distributions in order to get combined estimate
of age of surface"""

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
import matplotlib

#Rock Creek

data = np.genfromtxt('./data/test_ages.csv', delimiter=',', skip_header=1) [:,2:]
surface_labels = np.genfromtxt('./data/test_ages.csv', delimiter=',', skip_header=1, usecols=0, dtype=str)

# Get samples we want from each surface
# Edit the code below to add additional surfaces.
# Name must match surface labels in first column of data file
age_dict = {'S1':[]}#, 'T2':[], 'T3':[]}
sigma_dict = {'S1':[]}#, 'S2':[], 'S3':[]}
sample_dict = {'S1':[]}#, 'S2':[], 'S3':[]}
for i, label in enumerate(surface_labels):
    if label == 'S1':
        age_dict['S1'].append(data[i][0])
        sigma_dict['S1'].append(data[i][1])
#    elif label == 'S2':
#        age_dict['S2'].append(data[i][0])
#        sigma_dict['S2'].append(data[i][1]) 
#    elif label == 'S3':
#        age_dict['S3'].append(data[i][0])
#        sigma_dict['S3'].append(data[i][1])
    else:
        print('Sample from unknown surface')

# Check we've read it in correctly
print(age_dict)
print(sigma_dict)
for key, value in age_dict.items():
    for i, mu in enumerate(value):
        dnorm = norm(loc=mu, scale=sigma_dict[key][i])
        samples = dnorm.rvs(size=10000)
        sample_dict[key].append(samples)
    # Make an example plot
    if key=='S1':
        # Estimate bounds from data, may need to adjust manually (see line below)
        xvals = np.arange(min(value)-6*max(sigma_dict[key]), max(value)+6*max(sigma_dict[key]), 0.1)
#        xvals = np.arange(0, 40, 0.1)
        for i, mu in enumerate(value):
            dnorm = norm(loc=mu, scale=sigma_dict[key][i])
            plt.plot(xvals, dnorm.pdf(xvals), c='0.6')
            p1, = plt.fill(xvals, dnorm.pdf(xvals), facecolor='0.8', label='Individual sample') 
        mu_samples = np.mean(sample_dict[key])
        sig_samples = np.std(sample_dict[key])
        dnorm_samp = norm(loc=mu_samples, scale=sig_samples)
        plt.plot(xvals, dnorm_samp.pdf(xvals), c='0.3', zorder=11)
        p2, = plt.fill(xvals, dnorm_samp.pdf(xvals), facecolor='0.4',
                       zorder=10, label = 'Combined age')#, alpha=0.5)
        plt.xlabel('Age (ka)')
        plt.ylabel('Density')
        plt.legend(handles=[p1,p2])
        plt.savefig('plots/sampling_method_test_example.png', dpi=300)
        
print('\nTest data')
for key, value in sample_dict.items():
    mean = np.mean(value)
    sigma = np.std(value)
    print(key)
    print('Mean age (ka)', mean)
    print('Sigma (kyr)', sigma)
Sample age uncertainty distributions in order to get combined estimate of age of surface #Rock Creek # Get samples we want from each surface # Edit the code below to add additional surfaces. # Name must match surface labels in first column of data file #, 'T2':[], 'T3':[]} #, 'S2':[], 'S3':[]} #, 'S2':[], 'S3':[]} #    elif label == 'S2': #        age_dict['S2'].append(data[i][0]) #        sigma_dict['S2'].append(data[i][1]) #    elif label == 'S3': #        age_dict['S3'].append(data[i][0]) #        sigma_dict['S3'].append(data[i][1]) # Check we've read it in correctly # Make an example plot # Estimate bounds from data, may need to adjust manually (see line below) #        xvals = np.arange(0, 40, 0.1) #, alpha=0.5)