<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analysis of classification results
Objective
Step1:  Load original model
Step2:  Load sample classification results
The implemented classification method does not return a single best-fit model, but an ensemble of probable model (as it is an MCMC sampling from the posterior). As a first test, we will therefore import single models first and check the misclassification rate defined as
Step3:  Results of the classification do not necessarily contain the same ids as the units in the initial model. This seems to be the case here, as well. Re-sort
Step4:  Now remap results and compare again
Step6:  Combined analysis in a single function
Note
Step7:  Determine validity of uncertainty estimate
In addition to single model realisations, an esitmate of model uncertainty is calculated (this is, actually, also one of the main "selling points" of the paper). So, we will now check if the correct model is actually in the range of the estimated model uncertainty bounds (i.e.
Step8:  We now need to perform the remapping similar to before, but now for the probability fields
Step9:  Determination of misclassification statistics
Next step
Step10:  First, we again need to check the assignment of the units/ class ids
Step11:  We can now determine the misclassification for all results
<ASSISTANT_TASK:>
Python Code:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
import sys, os
import matplotlib.pyplot as plt
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
repo_path = os.path.realpath('../..')
import pynoddy.history
import numpy as np
%matplotlib inline
Explanation: Analysis of classification results
Objective: read back in the classification results and compare to original model
End of explanation
import pynoddy.output
reload(pynoddy.output)
output_name = "feature_out"
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('x', 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') # note: YlOrRd colourmap should be suitable for colorblindness!
Explanation: Load original model:
End of explanation
f_set1 = open("../../sandbox/jack/features_lowres-5 with class ID.csv").readlines()
f_set1[0]
# initialise classification results array
cf1 = np.empty_like(nout.block)
# iterate through results and append
for f in f_set1[1:]:
    fl = f.rstrip().split(",")
    cf1[int(fl[0]),int(fl[1]),int(fl[2])] = int(fl[6])
f_set1[2:6]
nout.plot_section('x', data = cf1,
                  colorbar = True, title="", layer_labels = range(5),
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') 
# compare to original model:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
nout.plot_section('x', ax = ax1, 
                  colorbar = False, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') # note: YlOrRd colourmap should be suitable for colorblindness!
nout.plot_section('x', data = cf1,ax = ax2,
                  colorbar = False, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') 
Explanation: Load sample classification results
The implemented classification method does not return a single best-fit model, but an ensemble of probable model (as it is an MCMC sampling from the posterior). As a first test, we will therefore import single models first and check the misclassification rate defined as:
$$\mbox{MCR} = \frac{\mbox{Number of misclassified voxels}}{\mbox{Total number of voxels}}$$
End of explanation
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1 = ax1.imshow(nout.block[15,:,:].transpose(), 
           interpolation = 'none', cmap = 'YlOrRd', origin = 'lower left')
plt.colorbar(im1)
im2 = ax2.imshow(cf1[15,:,:].transpose(), 
           interpolation = 'none', 
           cmap = 'YlOrRd', origin = 'lower left')
print np.unique(nout.block)
print np.unique(cf1)
# define id mapping from cluster results to original:
# id_mapping = {2:1, 3:2, 4:5, 5:3, 1:4}
# remapping for result 4:
# id_mapping = {4:5, 3:4, 1:3, 5:2, 2:1}
# remapping for result 5:
id_mapping = {2:5, 1:4, 3:3, 5:2, 4:1}
Explanation: Results of the classification do not necessarily contain the same ids as the units in the initial model. This seems to be the case here, as well. Re-sort:
End of explanation
def re_map(id_val):
    return id_mapping[id_val]
re_map_vect = np.vectorize(re_map)
cf1_remap = re_map_vect(cf1)
# compare to original model:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
nout.plot_section('x', ax = ax1, 
                  colorbar = False, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') # note: YlOrRd colourmap should be suitable for colorblindness!
nout.plot_section('x', data = cf1_remap, ax = ax2,
                  colorbar = False, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') 
feature_diff = (nout.block != cf1_remap)
nout.plot_section('x', data = feature_diff,
                  colorbar = False, title="Difference between real and matched model",
                  cmap = 'YlOrRd') 
# Calculate the misclassification:
np.sum(feature_diff) / float(nout.n_total)
# Export misclassification to VTK:
misclass = feature_diff.astype('int')
nout.export_to_vtk(vtk_filename = "misclass", data=misclass)
Explanation: Now remap results and compare again:
Note: create a vectorised function to enable a direct re-mapping of the entire array while keeping the structure!
End of explanation
def calc_misclassification(nout, filename):
    Calculate misclassification for classification results data stored in file
    
    **Arguments**:
        - *nout* = NoddyOutput: original model (Noddy object)
        - *filename* = filename (with path): file with classification results
    
    f_set1 = open(filename).readlines()
    # initialise classification results array
    cf1 = np.empty_like(nout.block)
    # iterate through results and append
    for f in f_set1[1:]:
        fl = f.rstrip().split(",")
        cf1[int(fl[0]),int(fl[1]),int(fl[2])] = int(fl[6])
    # remap ids
    cf1_remap = re_map_vect(cf1)
    # determine differences in class ids:
    feature_diff = (nout.block != cf1_remap)
    # Calculate the misclassification:
    misclass = np.sum(feature_diff) / float(nout.n_total)
    return misclass
filename = r"../../sandbox/jack/features_lowres-4 with class ID.csv"
calc_misclassification(nout, filename)
Explanation: Combined analysis in a single function
Note: function assumes correct EOL character in data file (check/ adjust with vi: %s/\r/\r/g)
Problem: remapping is unfortunatley not identical!
End of explanation
# f_set1 = open("../../sandbox/jack/features_lowres-6 with class ID and Prob.csv").readlines()
f_set1 = open("../../sandbox/jack/features_lowres-8 with Prob (weak Beta).csv").readlines()
f_set1[0]
# initialise classification results array
cf1 = np.empty_like(nout.block)
# Initialise probability array
probs = np.empty((5, cf1.shape[0], cf1.shape[1], cf1.shape[2]))
# iterate through results and append
for f in f_set1[1:]:
    fl = f.rstrip().split(",")
    i,j,k = int(fl[0]),int(fl[1]),int(fl[2])
    # cf1[i,j,k] = int(fl[6])
    for i2 in range(5):
        probs[i2,i,j,k] = float(fl[i2+6])
Explanation: Determine validity of uncertainty estimate
In addition to single model realisations, an esitmate of model uncertainty is calculated (this is, actually, also one of the main "selling points" of the paper). So, we will now check if the correct model is actually in the range of the estimated model uncertainty bounds (i.e.: if all voxets values from the original model actually have a non-zero probability in the estimated model)!
First step: load estimated class probabilities:
End of explanation
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1 = ax1.imshow(nout.block[15,:,:].transpose(), 
           interpolation = 'none', cmap = 'YlOrRd', origin = 'lower left')
plt.colorbar(im2)
im2 = ax2.imshow(probs[4,15,:,:].transpose(), 
           interpolation = 'none', 
           cmap = 'YlOrRd', origin = 'lower left')
# Note: map now ids from original model to probability fields in results:
prob_mapping = {4:0, 5:1, 3:2, 1:3, 2:4}
# Check membership for each class in original model
for i in range(1,6):
    tmp = np.ones_like(nout.block) * (nout.block==i)
    # test if voxels have non-zero probability by checking conjunction with zero-prob voxels
    prob_zero = probs[prob_mapping[i],:,:,:] == 0
    misidentified = np.sum(tmp * prob_zero)
    print i, misidentified
    
prob_zero = probs[prob_mapping[1],:,:,:] == 0
Explanation: We now need to perform the remapping similar to before, but now for the probability fields:
End of explanation
f_set1 = open("../../sandbox/jack/features_lowres-7 with 151 realizations.csv").readlines()
# Initialise results array
all_results = np.empty((152, cf1.shape[0], cf1.shape[1], cf1.shape[2]))
# iterate through results and append
for f in f_set1[1:]:
    fl = f.rstrip().split(",")
    i,j,k = int(fl[0]),int(fl[1]),int(fl[2])
    # cf1[i,j,k] = int(fl[6])
    for i2 in range(152):
        try:
            all_results[i2,i,j,k] = float(fl[i2+5])
        except IndexError:
            print i2, i, j, k
Explanation: Determination of misclassification statistics
Next step: use multiple results from one chain to determine misclassification statistics.
End of explanation
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1 = ax1.imshow(nout.block[15,:,:].transpose(), 
           interpolation = 'none', cmap = 'YlOrRd', origin = 'lower left')
plt.colorbar(im1)
im2 = ax2.imshow(all_results[5,15,:,:].transpose(), 
           interpolation = 'none', 
           cmap = 'YlOrRd', origin = 'lower left')
# mapping from results to original:
id_mapping = {2:5, 1:4, 3:3, 5:2, 4:1}
def re_map(id_val):
    return id_mapping[id_val]
re_map_vect = np.vectorize(re_map)
# Apply remapping to all but first result (seems to be original feature)
all_results_remap = re_map_vect(all_results[1:,:,:,:])
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1 = ax1.imshow(nout.block[30,:,:].transpose(), 
           interpolation = 'none', cmap = 'YlOrRd', origin = 'lower left')
# plt.colorbar(im1)
im2 = ax2.imshow(all_results_remap[85,30,:,:].transpose(), 
           interpolation = 'none', 
           cmap = 'YlOrRd', origin = 'lower left')
Explanation: First, we again need to check the assignment of the units/ class ids:
End of explanation
all_misclass = np.empty(151)
for i in range(151):
    # determine differences in class ids:
    feature_diff = (nout.block != all_results_remap[i,:,:,:])
    # Calculate the misclassification:
    all_misclass[i] = np.sum(feature_diff) / float(nout.n_total)
plt.plot(all_misclass)
Explanation: We can now determine the misclassification for all results:
End of explanation
<END_TASK>