<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 for drillholes in Gippsland Basin Model
We here evaluate how to analyse and visualise uncertainties in a kinematic model. The basic idea is that we have a set of drillhole locations and depths and want to know how uncertain the model is at these specific locations.
The required methods are implemented in an Experiment subclass and tested here with an application to the Gippsland Basin model.
Step1:  Creating an experiment object
First, we start with generating a pynoddy experiment object. The experiment class inherits all the methods from the base pynoddy.history class and we can directly import the Gippsland Basin model that we want to analyse into the object
Step2:  Some basic information about the model can be obtained with
Step3:  We can have a quick look at the model in a section view (note that Noddy is now executed in the background when required - and the output automatically generated in the required resolution)
Step4:  The base plot is not very useful - but we can create a section plot with a define vertical exaggeration (keyword ve) and plot the colorbar in horizontal orientation
Step5:  Note
Step6:  Generating random perturbations of the model
Before generating random prerturbations, we should now store the base version so that we can always revert to it at a later stage
Step7:  For a reproducible experiment, we can also set the random seed
Step8:  And now, let's perturb the model
Step9:  Let's see what happened
Step10:  ...and another perturbation
<ASSISTANT_TASK:>
Python Code:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
%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('/Users/flow/git/pynoddy/')
sys.path.append(repo_path)
import pynoddy.history
import pynoddy.experiment
rcParams.update({'font.size': 20})
Explanation: Uncertainty analysis for drillholes in Gippsland Basin Model
We here evaluate how to analyse and visualise uncertainties in a kinematic model. The basic idea is that we have a set of drillhole locations and depths and want to know how uncertain the model is at these specific locations.
The required methods are implemented in an Experiment subclass and tested here with an application to the Gippsland Basin model.
End of explanation
import importlib
importlib.reload(pynoddy.history)
importlib.reload(pynoddy.output)
importlib.reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file) 
Explanation: Creating an experiment object
First, we start with generating a pynoddy experiment object. The experiment class inherits all the methods from the base pynoddy.history class and we can directly import the Gippsland Basin model that we want to analyse into the object:
End of explanation
print(gipps_topo_ex)
Explanation: Some basic information about the model can be obtained with:
End of explanation
gipps_topo_ex.plot_section('y')
Explanation: We can have a quick look at the model in a section view (note that Noddy is now executed in the background when required - and the output automatically generated in the required resolution):
End of explanation
# gipps_topo_ex.determine_model_stratigraphy()
gipps_topo_ex.plot_section('x', ve = 5, position = 'centre',
                cmap = 'YlOrRd',
                title = '',
                colorbar = False)
gipps_topo_ex.plot_section('y', position = 100, ve = 5.,
                cmap = 'YlOrRd',
                title = '',
                colorbar_orientation = 'horizontal')
Explanation: The base plot is not very useful - but we can create a section plot with a define vertical exaggeration (keyword ve) and plot the colorbar in horizontal orientation:
End of explanation
importlib.reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file) 
gipps_topo_ex.load_parameter_file(os.path.join(repo_path, "examples/gipps_params_2.csv"))
Explanation: Note: The names of the model stratigraphy (colorbar labels) are unfortunately not defined correctly in the input file - we need to fix that, then we should get useful labels, as well!
Loading parameters ranges from file
We now need to define the parameter ranges. This step can either be done through explicit definition in the notebook (see the previous notebook on the Experiment class), or a list of parameters and defined ranges plus statistics can be read in from a csv file. This enables the convenient parameter definition in a spreadsheed (for example through Excel).
In order to be read in correctly, the header should contain the labels:
'event' : event id
'parameter' : Noddy parameter ('Dip', 'Dip Direction', etc.)
'min' : minimum value
'max' : maximum value
'initial' : initial value
In addition, it is possible to define PDF type and parameters. For now, the following settings are supported:
'type' = 'normal' 
'stdev' : standard deviation
'mean' : mean value (default: 'initial' value)
We can read in the parameters simply with:
End of explanation
gipps_topo_ex.freeze()
Explanation: Generating random perturbations of the model
Before generating random prerturbations, we should now store the base version so that we can always revert to it at a later stage:
End of explanation
gipps_topo_ex.set_random_seed(12345)
Explanation: For a reproducible experiment, we can also set the random seed:
End of explanation
gipps_topo_ex.random_perturbation()
Explanation: And now, let's perturb the model:
End of explanation
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base", 
                           colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False, 
                           title = "", ve = 5.)
#
# Note: keep these lines only for debugging!
#
reload(pynoddy.output)
reload(pynoddy.history)
reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file) 
gipps_topo_ex.load_parameter_file(os.path.join(repo_path, "examples/gipps_params.csv"))
# freeze base state
gipps_topo_ex.freeze()
# set seed
gipps_topo_ex.set_random_seed(12345)
# randomize
gipps_topo_ex.random_perturbation()
b1 = gipps_topo_ex.get_section('x', resolution = 50, model_type = 'base')
# b1.plot_section(direction = 'x', colorbar = False, title = "", ve = 5.)
b2 = gipps_topo_ex.get_section('x', resolution = 50, model_type = 'current')
# b1.plot_section(direction = 'x', colorbar = True, title = "", ve = 5.)
b1 -= b2
# b1.plot_section(direction = 'x', colorbar = True, title = "", ve = 5.)
print np.min(b1.block), np.max(b1.block)
type(b1)
Explanation: Let's see what happened: we can compare the new model to the base model as we stored it before:
End of explanation
gipps_topo_ex.random_perturbation()
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base", 
                           colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False, 
                           title = "", ve = 5.)
# plot difference
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base", 
                           colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False, 
                           title = "", ve = 5.)
gipps_topo_ex.param_stats
Explanation: ...and another perturbation:
End of explanation
<END_TASK>