<jupyter_start><jupyter_text>The workflow combines workflows in build_database_05.ipynb and prediction_03.ipynb workflows to test our latest model (CaCO3) on the core PS75-056-1. This core only has CaCO3 as the bulk chemistry.<jupyter_code>import numpy as np 
import pandas as pd
import glob
import matplotlib.pyplot as plt

#plt.style.use('ggplot')
plt.style.use('seaborn-colorblind')
plt.style.use('dark_background')
plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['savefig.bbox'] = 'tight'
plt.rcParams['savefig.transparent'] = True

%matplotlib inline

import datetime
date = datetime.datetime.now().strftime('%Y%m%d')<jupyter_output><empty_output><jupyter_text>Read and build spectral datasetThe core has the same data format as the usual one.<jupyter_code>file_name = []
spe_all = []
depth_all = []
cps_all = []
core_all = []
s_depth_all = []

# only read the 10kV which having better signal to the light elements
spe_dir = glob.glob('data/PS75-056-1/Run 1 at  10kV/*.spe')

# make sure the order follows the depthes in filename
spe_dir.sort()

for spe in spe_dir:
    check_depth = spe.split()[3].split('_')[-1]

    # there are some inconsistencies in nameing...as usual
    # 5 and 6 digis means in mm
    if len(check_depth) >= 5:
        start_depth = int(check_depth)
    # 3 and 4 digit means in cm, needs to be multipled to be mm
    elif len(check_depth) >= 3:
        start_depth = int(check_depth) * 10

    file_name.append(spe.split('/')[-1])

    with open(spe, 'r') as f:
        content = []
        lines = f.readlines()
        for line in lines[49:]:
            content = np.hstack((content, line.split()))

        section_depth = int(lines[13][:-3])
        spe_all.append(content.astype(int))
        cps_all.append(int(lines[28]))
        core_all.append('PS75-056-1')
        s_depth_all.append(section_depth)
        depth_all.append(section_depth + start_depth)
spe_df = pd.DataFrame(spe_all, columns = [str(_) for _ in range(2048)])
spe_df['cps'] = cps_all
spe_df['core'] = core_all
spe_df['composite_depth_mm'] = depth_all
spe_df['section_depth_mm'] = s_depth_all
spe_df['filename'] = file_name
spe_df
spe_df[spe_df.isnull().any(axis=1)]<jupyter_output><empty_output><jupyter_text>Build composite_id<jupyter_code>spe_df.composite_depth_mm.max()
composite_id = []
for core, depth in zip(spe_df.core, spe_df.composite_depth_mm):
    composite_id.append('{}_{:05}'.format(core, depth))
    
spe_df['composite_id'] = composite_id
spe_df = spe_df.set_index('composite_id')<jupyter_output><empty_output><jupyter_text>Drop duplicates<jupyter_code>len(spe_df.drop_duplicates('composite_id', keep = 'last'))<jupyter_output><empty_output><jupyter_text>There is no duplicates because the length is the same.  Build section<jupyter_code>section_all = []
# make sure the order follows the core and composite depth
spe_df.sort_values(by = 'composite_id', axis = 0, inplace = True)

# I assume every core scanned from section 0 so the first section in the core is marked as section 0
# the deeper the larger number
section = 0
X = spe_df['section_depth_mm']
for i in range(len(X)):
    section_all.append(section)
    try:
        # when section changes, the section depth should be rest to smaller number
        if X[i] > X[i + 1]:
            section += 1
    except IndexError:
        print('bottom of the core.')

spe_df['section'] = section_all
spe_df<jupyter_output><empty_output><jupyter_text>Read bulk chemistry<jupyter_code>bulk_df = pd.read_table('data/Bulk chem/PS75-56_bulkCaCO3 (%)', skiprows=1)
bulk_df['depth m'] = bulk_df['depth m'] * 1000
bulk_df.columns = ['mid_depth_mm', 'CaCO3%'] # modify headers to match previous bulk datasets
bulk_df['core'] = ['PS75-056-1' for _ in range(len(bulk_df))]
bulk_df = bulk_df[['CaCO3%', 'core', 'mid_depth_mm']]# modify orders to match previous bulk datasets
bulk_df
bulk_df.isna().any()<jupyter_output><empty_output><jupyter_text>Merge spe and bulk datasets<jupyter_code>mask_c = spe_df.columns[:2048]  # only the channels
merge_df = pd.DataFrame()

for index, row in bulk_df.iterrows():
    mid = row['mid_depth_mm']
    
    # get the spe in 10 mm interval
    mask_r = (spe_df.composite_depth_mm >= (mid-5)) & (spe_df.composite_depth_mm <= (mid+5))
    merge_df = pd.concat(
        [merge_df, spe_df.loc[mask_r, mask_c].apply(np.mean, axis = 0).append(row)],
        axis = 1
    )
    
merge_df = merge_df.T.reset_index(drop = True)
merge_df<jupyter_output><empty_output><jupyter_text>Check rows having nan in any column<jupyter_code>merge_df[merge_df.isnull().any(axis = 1)]<jupyter_output><empty_output><jupyter_text>Good, no nan in the dataset.  Export datasetThis dataset combines the merged datasets of core PS75-056-1 and only has CaCO3 measurements.<jupyter_code>merge_df.to_csv('data/spe+bulk_PS75-056-1_{}.csv'.format(date))
bulk_df.to_csv('data/bulk_PS75-056-1_{}.csv'.format(date))
spe_df.to_csv('data/spe_PS75-056-1_{}.csv'.format(date))
print(date)<jupyter_output>20210601<jupyter_text>Test the accuracy<jupyter_code>#merge_df = pd.read_csv('data/spe+bulk_PS75-056-1_20210601.csv', index_col=0)
X = merge_df.iloc[:, :2048].values
X = X / X.sum(axis = 1, keepdims = True)
from joblib import load
m_caco3 = load('models/caco3_nmf+svr_model_20201216.joblib')<jupyter_output><empty_output><jupyter_text>Predict<jupyter_code>from sklearn.metrics import r2_score
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import max_error

y_caco3 = np.exp(m_caco3.predict(X))
print('Scores in the test set:')
print('R2 = {:.3f} .'.format(r2_score(merge_df['CaCO3%'], y_caco3)))
print('The mean absolute error is {:.3f} (%, concetration).'.format(mean_absolute_error(merge_df['CaCO3%'], y_caco3)))
print('The max. residual error is {:.3f} (%, concetration).'.format(max_error(merge_df['CaCO3%'], y_caco3)))
plt.plot(range(len(y_caco3)), merge_df['CaCO3%'], alpha=0.6, label='Measurement')
plt.plot(range(len(y_caco3)), y_caco3, label='Prediction (R$^2$={:.2f})'.format(r2_score(merge_df['CaCO3%'], y_caco3)))
#plt.text(12, -7, r'R$^2$={:.2f}, mean ab. error={:.2f}, max. ab. error={:.2f}'.format(grid.best_score_, mean_absolute_error(y_ttest, y_predict), max_error(y_ttest, y_predict)))
plt.ylabel('CaCO$_3$ concentration (%)')
plt.xlabel('Sample no.')
plt.legend(loc = 'upper right')
plt.savefig('results/caco3_predictions_PS75-056-1_{}.png'.format(date))
plt.hist(merge_df['CaCO3%'])
plt.xlabel('CaCO${_3}$ (wt%) ')
plt.ylabel('Count')
plt.savefig('results/caco3_hist_PS75-056-1_{}.png'.format(date))
plt.plot(range(520), X.max(axis=0)[:520], label='Max.')
plt.plot(range(520), X.min(axis=0)[:520], label='Min.')
fig, ax = plt.subplots(1, 1)
ax.plot([0, 80], [0, 80], ls='--', c='gray')
ax.scatter(merge_df['CaCO3%'], y_caco3, s=6)
ax.set_aspect('equal')
ax.set_xlabel('Measurement (wt%) ')
ax.set_ylabel('Prediction (wt%) ')
fig.savefig('results/caco3_predictions_PS75-056-1_2_{}.png'.format(date))<jupyter_output><empty_output><jupyter_text>Draw spectrum<jupyter_code>pre_df = pd.read_csv('data/spe+bulk_dataset_20201215.csv', index_col=0)
pre_df.core.unique()
from sklearn.model_selection import train_test_split
X_p = pre_df.iloc[:, : -5].values
X_p = X_p / X_p.sum(axis = 1, keepdims = True)

SO264_core = pre_df.core.unique()[:-3]

s1_X_train, _ = train_test_split(X_p[pre_df.core.isin(SO264_core)], test_size = 0.2, shuffle = True, random_state = 24)
s2_ca_train, _ = train_test_split(X_p[pre_df.core.isin(['LV28-44-3', 'LV29-114-3'])], test_size = 0.2, shuffle = True, random_state = 24)<jupyter_output><empty_output><jupyter_text>The sum of carbonate and TOC training data amounts fit to the histogram, respectively. I don't have enough time to fine-tune a figure for all these spectra in different steps and analytes so I plot the figure of the 1st+2nd scale ups only for now.<jupyter_code>fig, axes = plt.subplots(1, 2, figsize=(7.2, 3), sharex='row', sharey='row')

axes[0].plot(range(2048), np.vstack((s1_X_train, s2_ca_train)).max(axis=0), label='Max.')
axes[0].plot(range(2048), np.vstack((s1_X_train, s2_ca_train)).min(axis=0), label='Min.')
axes[0].set_xlim(0, 520)
axes[0].set_ylabel('Normalized intensity')
axes[0].set_xlabel('Channel, Training data')

axes[1].plot(range(2048), X.max(axis=0), label='Max.')
axes[1].plot(range(2048), X.min(axis=0), label='Min.')
axes[1].set_xlabel('Channel, PS75-056-1')
axes[1].legend()

fig.subplots_adjust(wspace=.04)
plt.savefig('results/spectra_comparison_1_{}.png'.format(date))
print(date)
plt.plot(range(2048), np.vstack((s1_X_train, s2_ca_train)).mean(axis=0), label='Training data')
plt.plot(range(2048), X.mean(axis=0), label='PS75-056-1')
plt.xlim(0, 520)
plt.ylabel('Normalized intensity')
plt.xlabel('Channel')
plt.legend()
plt.suptitle('Averaged spectra')
plt.subplots_adjust(top=.92)
plt.savefig('results/spectra_comparison_2_{}.png'.format(date))
print(date)<jupyter_output>20210603