<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Uncertainty analysis of a 2-D slice model
Possible paper titles
Step1:  Model set-up
Subsequently, we will use a model from the "Atlas of Structural Geophysics" as an example model.
Step2:  BUG!!!!
Note
Step3:  We now define the parameter uncertainties
Step4:  And, in a next step, perform the model sampling
Step5:  Save the model data for later re-use (e.g. to extend the data set)
Step7:  Testing parallel execution
As a next step, use parallel execution of noddy calculation - should be relatively simple, however
Step8:  Adapt model generation to use temp directory
Step9:  Try compute_iter_pool function in normal framework (non-parallel execution)
Step10:  Trying with pool
Step11:  Executing Noddy MC in parallel
After a bit of fiddling around, it seems to work now
Step12:  Use hspace to calculate joint entropies
Adjusted notebook
Step13:  Calculation of cell information entropy
(Include note on
Step14:  We now visualise the cell information entropy, shown in Fig. (). We can here clearly identify uncertain regions within this model section. It is interesting to note that we can mostly still identify the distinct layer boundaries in the fuzzy areas of uncertainty around their borders (note
Step15:  Here again an example of single models (adjust to visualise and probably include something like a "table plot" of multiple images for a paper!)
Step16:  And here the "mean" lithologies (note
Step17:  And here a bit more meaningful
Step18:  Idea
Step19:  Try own "entropy colormap"
Step20:  For comparison again
Step21:  And the difference, for clarity
Step22:  Clearly, the highset reduction is in the area around the borehole, but interestingly, the uncertianty in other areas is also reduced! Note specifically the reduction of uncertainties in the two neighbouring fold hinges.
Let's check some other positions (and drilling "depths")
Step23:  We also just include one timing step to estimate the approximate simualtion time
Step24:  Intersting! Only a local reduction around the drilling position, however
Step25:  Check deep "drilling" at pos 60
Step26:  Interesting! And now both combined
Step27:  We can see that now only a part on the left remains with significant uncertainties. So, let's "drill" into this, as well
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
# here the usual imports. If any of the imports fails, 
# make sure that pynoddy is installed
# properly, ideally with 'python setup.py develop' 
# or 'python setup.py install'
import sys, os
import matplotlib.pyplot as plt
import numpy as np
# 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('../..')
sys.path.append('../..')
import pynoddy
import importlib
importlib.reload(pynoddy)
import pynoddy.history
import pynoddy.experiment
importlib.reload(pynoddy.experiment)
rcParams.update({'font.size': 15})
pynoddy.history.NoddyHistory(history="typeb.his")
Explanation: Uncertainty analysis of a 2-D slice model
Possible paper titles: 
Posterior analysis of geological models/ geophysical inversions? with information theory
...with measures from information theory
...with information theoretic measures
...with information measures
...with information
Ensemble analysis...
Posterior and ensemble analysis...
(Note: reference to posterior analysis, e.g. in Tarantola paper!)
Include:
Analysis of pynoddy models:
simple example of graben (because the uncertainty reduction is so counter-intuitive), 
and maybe more complex example of folding structure?
Analysis of object modelling results? (i.e. the typical "Stanford channels"? but: only two outcomes, so not so meaningful... better example?)
Analysis of posterior ensemble from geophysical inversion (Greenstone model? Other examples from Mark & Mark?)
Journal? Math. Geo? Tectonophysics (note: relevance to strucural geological models!)
Include theory: error bounds on information measures!
End of explanation
# from pynoddy.experiment import monte_carlo
model_url = 'http://tectonique.net/asg/ch3/ch3_7/his/typeb.his'
ue = pynoddy.experiment.Experiment(history="typeb.his")
import sys
if sys.version_info[0] < 3:
    raise "Must be using Python 3"
print(sys.version_info[1])
ue.change_cube_size(100)
sec = ue.get_section('y')
tmp = open("typeb.his").readlines()
Explanation: Model set-up
Subsequently, we will use a model from the "Atlas of Structural Geophysics" as an example model.
End of explanation
sec.block.shape
ue.plot_section('y')
plt.imshow(sec.block[:,50,:].transpose(), origin = 'lower left', interpolation = 'none')
tmp = sec.block[:,50,:]
tmp.shape
ue.set_random_seed(12345)
ue.info(events_only = True)
Explanation: BUG!!!!
Note: there is either a bug in pynoddy or in Noddy itself: but the slice plotting method fails: actually, not a slice is computed but the entire model (and the extent is also not correct!). Check with Mark in course of paper prep!
End of explanation
param_stats = [{'event' : 2, 
              'parameter': 'Amplitude',
              'stdev': 100.0,
              'type': 'normal'},
              {'event' : 2, 
              'parameter': 'Wavelength',
              'stdev': 500.0,
              'type': 'normal'},
              {'event' : 2, 
              'parameter': 'X',
              'stdev': 500.0,
              'type': 'normal'}]
ue.set_parameter_statistics(param_stats)
Explanation: We now define the parameter uncertainties:
End of explanation
ue.set_random_seed(112358)
# perfrom random sampling
resolution = 100
sec = ue.get_section('y')
n_draws = 100
tmp = sec.block[:,50,:]
#
# Note: setting the dtype to 'int8' significantly reduces file size!
#
model_sections = np.empty((n_draws, tmp.shape[0], tmp.shape[1]), dtype='int8')
for i in range(n_draws):
    ue.random_draw()
    tmp_sec = ue.get_section('y', resolution = resolution, 
                             remove_tmp_files = True)
    model_sections[i,:,:] = tmp_sec.block[:,50,:]
Explanation: And, in a next step, perform the model sampling:
End of explanation
import pickle
model_sections.dtype
import pickle
f_out = open("model_sections_tmp_int.pkl", 'wb')
pickle.dump(model_sections, f_out)
f_out.close()
f_in = open("model_sections_tmp_int.pkl", 'rb')
model_sections = pickle.load(f_in)
Explanation: Save the model data for later re-use (e.g. to extend the data set):
End of explanation
#
# Store current directory to get back from temporary files
#
ori_dir = os.getcwd()
import multiprocessing as mp
import random
import string
random.seed(123)
# Define an output queue
output = mp.Queue()
# define a example function
def rand_string(length, output):
     Generates a random string of numbers, lower- and uppercase chars. 
    rand_str = ''.join(random.choice(
                        string.ascii_lowercase
                        + string.ascii_uppercase
                        + string.digits)
                   for i in range(length))
    output.put(rand_str)
# Setup a list of processes that we want to run
processes = [mp.Process(target=rand_string, args=(5, output)) for x in range(4)]
# Run processes
for p in processes:
    p.start()
# Exit the completed processes
for p in processes:
    p.join()
# Get process results from the output queue
results = [output.get() for p in processes]
print(results)
Explanation: Testing parallel execution
As a next step, use parallel execution of noddy calculation - should be relatively simple, however: note that, potentially, tmp-files may be overwritten! This may actually require an adaptation of compute_model to run in a tmp-folder...
Anyway, as a first step: here an example for parallel execution adapted from
http://sebastianraschka.com/Articles/2014_multiprocessing.html#the-process-class
End of explanation
import tempfile
import shutil
# Define function to perform one iteration
# Execute iterations in temporary directories to avoid overlap
#
# Note: needs to take outout as argument to add results
def compute_iter(ue, output):
    ue.random_draw()
    dirpath = tempfile.mkdtemp()
    os.chdir(dirpath)
    tmp_sec = ue.get_section('y', resolution = 100, 
                             remove_tmp_files = True)
    output.put(tmp_sec.block[:,50,:])
    
# Note: this is not the case for the 'pool' method:
# use `with` context management method to ensure that directory is deleted afterwards:
def compute_iter_pool(ue):
    from tempfile import TemporaryDirectory
    with TemporaryDirectory() as temp_dir:
        os.chdir(temp_dir)
        np.random.seed()
        ue.random_draw()
        tmp_sec = ue.get_section('y', resolution = 100, 
                                 remove_tmp_files = True)
    
    return tmp_sec.block[:,50,:]
np.random.seed()
np.random.randint(10)
Explanation: Adapt model generation to use temp directory:
End of explanation
ori_dir = "/Users/Shared/git/pynoddy/docs/notebooks/"
os.chdir("/Users/Shared/git/pynoddy/docs/notebooks/") 
ue.set_random_seed(112358)
# perfrom random sampling
resolution = 100
sec = ue.get_section('y')
tmp = sec.block[:,50,:]
n_draws = 100
#
# Note: setting the dtype to 'int8' significantly reduces file size!
#
model_sections = np.empty((n_draws, tmp.shape[0], tmp.shape[1]), dtype='int8')
#
#
for i in range(n_draws):
    model_sections[i,:,:] = compute_iter_pool(ue)
os.chdir("/Users/Shared/git/pynoddy/docs/notebooks/")   
model_sections.shape
compute_iter_pool(ue)
os.chdir(ori_dir)
!pwd
Explanation: Try compute_iter_pool function in normal framework (non-parallel execution):
End of explanation
def cube(x):
    return x**3
pool = mp.Pool(processes=4)
results = [pool.apply(cube, args=(x,)) for x in range(1,7)]
print(results)
pool = mp.Pool(processes=4)
results = [pool.apply_async(cube, args=(ue,)) for x in range(1,7)]
output = [p.get() for p in results]
print(output)
Explanation: Trying with pool
End of explanation
import copy
ue.set_random_seed(112358)
# perfrom random sampling
resolution = 100
sec = ue.get_section('y')
tmp = sec.block[:,50,:]
n_draws = 100000
#
# Note: setting the dtype to 'int8' significantly reduces file size!
#
model_sections = np.empty((n_draws, tmp.shape[0], tmp.shape[1]), dtype='int8')
result_list = []
def log_result(result):
    # This is called whenever foo_pool(i) returns a result.
    # result_list is modified only by the main process, not the pool workers.
    result_list.append(result)
pool = mp.Pool(processes=20)
for i in range(n_draws):
    # pool.apply_async(foo_pool, args = (i, ), callback = log_result)
    # try to use copy, but this did not fix the problem, unfortunately...
    # ue_copy = copy.deepcopy(ue)
    pool.apply_async(compute_iter_pool, args=(ue,), callback = log_result)
pool.close()
pool.join()
# model_sections = np.array([pool.apply(compute_iter_pool, args=(ue,)) for x in range(4)])
ori_dir = "/Users/Shared/git/pynoddy/docs/notebooks/"
os.chdir(ori_dir)
model_sections = np.array(result_list)
model_sections.shape
import pickle
f_out = open("model_sections_tmp_int.pkl", 'wb')
pickle.dump(model_sections, f_out)
f_out.close()
model_sections[0,:,:]
plt.imshow(model_sections[1,:,:].transpose()-model_sections[0,:,:].transpose(), origin = 'lower left', interpolation = 'none')
ori_dir = os.getcwd()
print(ori_dir)
Explanation: Executing Noddy MC in parallel
After a bit of fiddling around, it seems to work now:
End of explanation
sys.path.append('/Users/Shared/git/hspace')
import hspace.measures
importlib.reload(hspace.measures)
model_sections.shape
Explanation: Use hspace to calculate joint entropies
Adjusted notebook: use hspace-package to calculate information theoretic measures
End of explanation
hspace.measures.joint_entropy(model_sections[:,50,:], [0,1])
h = np.empty_like(model_sections[0,:,:], dtype='float64')
for i in range(100):
    for j in range(40):
        h[i,j] = hspace.measures.joint_entropy(model_sections[:,i,j])
h[50,30]
Explanation: Calculation of cell information entropy
(Include note on: theory of entropy calculation)
(Include in this paper: estimates on error bounds?)
Here now the function to calculate entropy from a data array in general. What we will need to do later is to pass all results at a single position as a "data array" and we can then estimate the information entropy at this position.
This function already expects a sorted array as an input and then uses the (ultra-fast) switchpoint method to calculate entropy:
The algorithm works on the simple idea that we do not explicitly require the single outputs at each location, but only the relative probability values. This may not matter too much for single entropy estimates (uni-variate), but it will matter a lot for multivariate cases, because we do not need to check all possible outcomes! Note that all outcomes with zero probability are simply not considered in the sorting algorithm (and they do not play any role in the calculation of the entropy, anyway), and that's exactly what we want to have!
In this new version, we use the implementation the hspace package:
End of explanation
plt.imshow(h.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
Explanation: We now visualise the cell information entropy, shown in Fig. (). We can here clearly identify uncertain regions within this model section. It is interesting to note that we can mostly still identify the distinct layer boundaries in the fuzzy areas of uncertainty around their borders (note: highlight in Figure!). However, additional aspects of uncertainty are now introduced: (a) the uncertainty about the x-position of the folds (see parameters: event 2, parameter x) is now clearly visible, and (b) uncertianties now seem to concentrate on the fold hinges. However, this is not so clear in the left part of the model, where the fold hing seems to be the least uncertain part. (check why: is this where the fold is actually fixed (even though still uncertain). My current interpretation: the fold location is fixed somewhere near this point, and so the wavelength uncertainty does not play a significant role. Furthermore, the fold is quite "open" at this position (i.e.: low angle between hinges) and therefore lateral shifts do not play a significant role.
End of explanation
plt.imshow(model_sections[1,:,:].transpose(), origin = 'lower left', interpolation = 'none')
Explanation: Here again an example of single models (adjust to visualise and probably include something like a "table plot" of multiple images for a paper!):
End of explanation
plt.imshow(np.mean(model_sections, axis = 0).transpose(), origin = 'lower left', interpolation = 'none')
Explanation: And here the "mean" lithologies (note: not really a distinct meaning, simply showing the average litho ids - could be somehow interpreted as characteristic functions, though...).
End of explanation
# step 1: estimate probabilities (note: unfortunate workaround with ones multiplication,
# there may be a better way, but this is somehow a recurring problem of implicit
# array flattening in numpy)
litho_id = 2
prob = np.sum(np.ones_like(model_sections) * (model_sections == litho_id), axis = 0) / model_sections.shape[0]
plt.imshow(prob.transpose(), 
           origin = 'lower left', 
           interpolation = 'none', 
           cmap = 'gray_r')
plt.colorbar(orientation = 'horizontal')
Explanation: And here a bit more meaningful: the analysis of single layer probabilities:
End of explanation
importlib.reload(hspace.measures)
dx = 15
xvals = np.ones(dx, dtype=int) * 20
yvals = np.arange(11,11+dx, dtype=int)
pos = np.vstack([xvals, yvals])
hspace.measures.joint_entropy(model_sections, pos.T)
# now: define position of "drill":
nx = 10
xvals = np.ones(nx, dtype=int) * 60
yvals = np.arange(39,39-nx, -1, dtype=int)
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add new position to positions vector:
        pos_all = np.vstack([pos, np.array([i,j])])
        # determine joint entropy
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow(h_cond_drill.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,100])
plt.ylim([0,40])
Explanation: Idea: also include a "new" consideration: where to collect information to reduce uncertainty of a single layer? Could be identified by reducing layer fuzziness, for example!
Or: what are most likely positions/ locations of a specific unit, given collected information?
<div class="alert alert-warning">
<b>More ideas:</b> 
<ul>
<li>General methods to create simple representations in 1-D, 2-D?
<li>Automatic probability plots (note: simple in 1-D, but also for slices in 2-D?) => Compare visualisations in previous notebooks!
</ul>
</div>
Analysis of multivariate condtional entropy
Later also:
"opposite" question, i.e.: if we would like to resolve uncertainty in a specific region: where to look best?
"complete" uncertainty (i.e.: joint entropy!)
greedy search for best spots for uncertainty reduction, simple (cell-wise), complex (related to potential drilling positions)
further ideas for "greedy search" to reduce uncertainties in a specific "search region" (i.e.: expected location of a deposit, etc.):
start with cell with highest (multivariate) mutual information
rank cells with highest I due to their own mutual information with other cells, which are not part of a defined "search region"
for simple, cell-wise: describe similarity to mapping! Maybe even a field example with data from Belgium? But this is still one or two MSc theses away...)
For the joint entropy analysis, we now use the new lexicographic (correct term?) sorting algorithm, implemented in the module hspace:
End of explanation
n_max = 5
h_max = np.max(h)
plt.imshow(h_cond_drill.transpose(), origin = 'lower left',
          cmap = 'viridis', interpolation = 'none', vmax=h_max)
plt.colorbar(orientation = 'horizontal')
# half-step contour lines
contour_levels = np.log2(np.arange(1., n_max + 0.001, .5))
plt.contour(h_cond_drill.transpose(), contour_levels, colors = 'gray')
# superpose 1-step contour lines
contour_levels = np.log2(np.arange(1., n_max + 0.001, 1.))
plt.contour(h_cond_drill.transpose(), contour_levels, colors = 'white')
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: Try own "entropy colormap":
End of explanation
plt.imshow(h.transpose(), origin = 'lower left',
          cmap = 'viridis', interpolation = 'none', vmax=h_max)
plt.colorbar(orientation = 'horizontal')
# half-step contour lines
contour_levels = np.log2(np.arange(1., n_max + 0.001, .5))
plt.contour(h.transpose(), contour_levels, colors = 'gray')
# superpose 1-step contour lines
contour_levels = np.log2(np.arange(1., n_max + 0.001, 1.))
plt.contour(h.transpose(), contour_levels, colors = 'white')
Explanation: For comparison again: the entropy of the initial model:
End of explanation
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'viridis', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: And the difference, for clarity:
End of explanation
# define position of "drill":
nx = 10
xvals = np.ones(nx, dtype=int) * 20
yvals = np.arange(39,39-nx, -1, dtype=int)
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
Explanation: Clearly, the highset reduction is in the area around the borehole, but interestingly, the uncertianty in other areas is also reduced! Note specifically the reduction of uncertainties in the two neighbouring fold hinges.
Let's check some other positions (and drilling "depths"):
End of explanation
import time
start = time.time()
pos_all = np.vstack([pos, np.array([50,20])])
h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
end = time.time()
# print(end-start)
# esimated total time:
ttime = 100 * 40 * (end-start)
print("Estimated total time: %.3f seconds or %.3f minutes" % (ttime, ttime/60.))
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add position to locations
        pos_all = np.vstack([pos, np.array([i,j])])
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: We also just include one timing step to estimate the approximate simualtion time:
End of explanation
# define position of "drill":
nx = 30
xvals = np.ones(nx, dtype=int) * 20
yvals = np.arange(39,39-nx, -1, dtype=int)
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add position to locations
        pos_all = np.vstack([pos, np.array([i,j])])
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
plt.imshow(h.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: Intersting! Only a local reduction around the drilling position, however: extending to the deeper layers, as well! Why?
Drill deeper:
End of explanation
# define position of "drill":
nx = 30
xvals = np.ones(nx, dtype=int) * 60
yvals = np.arange(39,39-nx, -1, dtype=int)
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add position to locations
        pos_all = np.vstack([pos, np.array([i,j])])
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow(h_cond_drill.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: Check deep "drilling" at pos 60
End of explanation
# define position of "drill":
nx = 30
xvals = np.hstack([np.ones(nx, dtype=int) * 60, np.ones(nx, dtype=int) * 30])
yvals = np.hstack([np.arange(39,39-nx, -1, dtype=int), np.arange(39,39-nx, -1, dtype=int)])
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
%%timeit
h_joint_loc = hspace.measures.joint_entropy(model_sections, pos)
# esimated total time:
ttime = 100 * 40 * 0.0002
print("Estimated total time: %.3f seconds or %.3f minutes" % (ttime, ttime/60.))
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add position to locations
        pos_all = np.vstack([pos, np.array([i,j])])
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow(h_cond_drill.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: Interesting! And now both combined:
End of explanation
# define position of "drill":
nx = 30
xvals = np.hstack([np.ones(nx, dtype=int) * 60, 
                   np.ones(nx, dtype=int) * 30,
                   np.ones(nx, dtype=int) * 5])
yvals = np.hstack([np.arange(39,39-nx, -1, dtype=int), 
                   np.arange(39,39-nx, -1, dtype=int), 
                   np.arange(39,39-nx, -1, dtype=int)])
pos = np.vstack([xvals, yvals]).T
# determine joint entropy of drill_locs:
h_joint_drill = hspace.measures.joint_entropy(model_sections, pos)
# generate conditional entropies for entire section:
h_cond_drill = np.zeros_like(h)
for i in range(100):
    for j in range(40):
        # add position to locations
        pos_all = np.vstack([pos, np.array([i,j])])
        h_joint_loc = hspace.measures.joint_entropy(model_sections, pos_all)
        # subtract joint entropy of drill locs to obtain conditional entropy
        h_cond_drill[i,j] = h_joint_loc - h_joint_drill
plt.imshow(h_cond_drill.transpose(), origin = 'lower left',
          cmap = 'gray', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
plt.imshow((h - h_cond_drill).transpose(), origin = 'lower left',
          cmap = 'RdBu', interpolation = 'none')
plt.colorbar(orientation = 'horizontal')
# plot drilling positions above it:
plt.plot(pos[:,0], pos[:,1], 'ws')
plt.xlim([0,99])
plt.ylim([0,39])
Explanation: We can see that now only a part on the left remains with significant uncertainties. So, let's "drill" into this, as well:
End of explanation
<END_TASK>