docs

docs/_downloads/066a96b3e7333f9a789eb7e687c45875/plot_histogram_2d.py

@@ -0,0 +1,305 @@
+"""
+Histogram 2D
+------------
+
+This 2D histogram class allows efficient updating of histograms, plotting and
+saving as HDF5.
+
+"""
+
+#%%
+import h5py
+import geobipy
+from geobipy import StatArray
+from geobipy import Histogram
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from geobipy import RectilinearMesh2D
+import numpy as np
+
+
+#%%
+# Create some histogram bins in x and y
+x = StatArray(np.linspace(-4.0, 4.0, 100), 'Variable 1')
+y = StatArray(np.linspace(-4.0, 4.0, 105), 'Variable 2')
+
+mesh = RectilinearMesh2D(x_edges=x, y_edges=y)
+#%%
+# Instantiate
+H = Histogram(mesh)
+
+#%%
+# Generate some random numbers
+a = np.random.randn(1000000)
+b = np.random.randn(1000000)
+
+#%%
+# Update the histogram counts
+H.update(a, b)
+
+#%%
+plt.figure()
+plt.subplot(131)
+plt.title("2D Histogram")
+_ = H.plot(cmap='gray_r')
+plt.subplot(132)
+H.pdf.plot(cmap='gray_r')
+plt.subplot(133)
+H.pmf.plot(cmap='gray_r')
+
+
+plt.figure()
+plt.subplot(131)
+H.cdf(axis=0).plot()
+plt.subplot(132)
+H.cdf(axis=1).plot()
+plt.subplot(133)
+H.cdf().plot()
+
+#%%
+# We can overlay the histogram with its credible intervals
+plt.figure()
+plt.title("90% credible intervals overlain")
+H.pcolor(cmap='gray_r')
+H.plotCredibleIntervals(axis=0, percent=95.0)
+_ = H.plotCredibleIntervals(axis=1, percent=95.0)
+
+#%%
+# Generate marginal histograms along an axis
+h1 = H.marginalize(axis=0)
+h2 = H.marginalize(axis=1)
+
+#%%
+# Note that the names of the variables are automatically displayed
+plt.figure()
+plt.suptitle("Marginals along each axis")
+plt.subplot(121)
+h1.plot()
+plt.subplot(122)
+_ = h2.plot()
+
+#%%
+# Create a combination plot with marginal histograms.
+# sphinx_gallery_thumbnail_number = 3
+plt.figure()
+gs = gridspec.GridSpec(5, 5)
+gs.update(wspace=0.3, hspace=0.3)
+ax = [plt.subplot(gs[1:, :4])]
+H.pcolor(colorbar = False)
+
+ax.append(plt.subplot(gs[:1, :4]))
+h = H.marginalize(axis=0).plot()
+plt.xlabel(''); plt.ylabel('')
+plt.xticks([]); plt.yticks([])
+ax[-1].spines["left"].set_visible(False)
+
+ax.append(plt.subplot(gs[1:, 4:]))
+h = H.marginalize(axis=1).plot(transpose=True)
+plt.ylabel(''); plt.xlabel('')
+plt.yticks([]); plt.xticks([])
+ax[-1].spines["bottom"].set_visible(False)
+
+#%%
+# Take the mean or median estimates from the histogram
+mean = H.mean()
+median = H.median()
+
+#%%
+plt.figure(figsize=(9.5, 5))
+plt.suptitle("Mean, median, and credible interval overlain")
+ax = plt.subplot(121)
+H.pcolor(cmap='gray_r', colorbar=False)
+H.plotCredibleIntervals(axis=0)
+H.plotMedian(axis=0, color='g')
+H.plotMean(axis=0, color='y')
+plt.legend()
+
+plt.subplot(122, sharex=ax, sharey=ax)
+H.pcolor(cmap='gray_r', colorbar=False)
+H.plotCredibleIntervals(axis=1)
+H.plotMedian(axis=1, color='g')
+H.plotMean(axis=1, color='y')
+plt.legend()
+
+#%%
+# Get the range between credible intervals
+H.credible_range(percent=95.0)
+
+#%%
+# We can map the credible range to an opacity or transparency
+H.opacity()
+H.transparency()
+
+# H.animate(0, 'test.mp4')
+
+import h5py
+with h5py.File('h2d.h5', 'w') as f:
+    H.toHdf(f, 'h2d')
+
+with h5py.File('h2d.h5', 'r') as f:
+    H1 = Histogram.fromHdf(f['h2d'])
+
+plt.close('all')
+
+x = StatArray(5.0 + np.linspace(-4.0, 4.0, 100), 'Variable 1')
+y = StatArray(10.0 + np.linspace(-4.0, 4.0, 105), 'Variable 2')
+
+mesh = RectilinearMesh2D(x_edges=x, x_relative_to=5.0, y_edges=y, y_relative_to=10.0)
+#%%
+# Instantiate
+H = Histogram(mesh)
+
+#%%
+# Generate some random numbers
+a = np.random.randn(1000000) + 5.0
+b = np.random.randn(1000000) + 10.0
+
+#%%
+# Update the histogram counts
+H.update(a, b)
+
+#%%
+plt.figure()
+plt.subplot(131)
+plt.title("2D Histogram")
+_ = H.plot(cmap='gray_r')
+plt.subplot(132)
+H.pdf.plot(cmap='gray_r')
+plt.subplot(133)
+H.pmf.plot(cmap='gray_r')
+
+plt.figure()
+plt.subplot(131)
+H.cdf(axis=0).plot()
+plt.subplot(132)
+H.cdf(axis=1).plot()
+plt.subplot(133)
+H.cdf().plot()
+
+#%%
+# We can overlay the histogram with its credible intervals
+plt.figure()
+plt.title("90% credible intervals overlain")
+H.pcolor(cmap='gray_r')
+H.plotCredibleIntervals(axis=0, percent=95.0)
+_ = H.plotCredibleIntervals(axis=1, percent=95.0)
+
+# Generate marginal histograms along an axis
+h1 = H.marginalize(axis=0)
+h2 = H.marginalize(axis=1)
+
+#%%
+# Note that the names of the variables are automatically displayed
+plt.figure()
+plt.suptitle("Marginals along each axis")
+plt.subplot(121)
+h1.plot()
+plt.subplot(122)
+_ = h2.plot()
+
+#%%
+# Create a combination plot with marginal histograms.
+# sphinx_gallery_thumbnail_number = 3
+plt.figure()
+gs = gridspec.GridSpec(5, 5)
+gs.update(wspace=0.3, hspace=0.3)
+ax = [plt.subplot(gs[1:, :4])]
+H.pcolor(colorbar = False)
+
+ax.append(plt.subplot(gs[:1, :4]))
+h = H.marginalize(axis=0).plot()
+plt.xlabel(''); plt.ylabel('')
+plt.xticks([]); plt.yticks([])
+ax[-1].spines["left"].set_visible(False)
+
+ax.append(plt.subplot(gs[1:, 4:]))
+h = H.marginalize(axis=1).plot(transpose=True)
+plt.ylabel(''); plt.xlabel('')
+plt.yticks([]); plt.xticks([])
+ax[-1].spines["bottom"].set_visible(False)
+
+#%%
+# Take the mean or median estimates from the histogram
+mean = H.mean()
+median = H.median()
+
+#%%
+plt.figure(figsize=(9.5, 5))
+plt.suptitle("Mean, median, and credible interval overlain")
+ax = plt.subplot(121)
+H.pcolor(cmap='gray_r', colorbar=False)
+H.plotCredibleIntervals(axis=0)
+H.plotMedian(axis=0, color='g')
+H.plotMean(axis=0, color='y')
+plt.legend()
+
+plt.subplot(122, sharex=ax, sharey=ax)
+H.pcolor(cmap='gray_r', colorbar=False)
+H.plotCredibleIntervals(axis=1)
+H.plotMedian(axis=1, color='g')
+H.plotMean(axis=1, color='y')
+plt.legend()
+
+#%%
+# Get the range between credible intervals
+H.credible_range(percent=95.0)
+
+#%%
+# We can map the credible range to an opacity or transparency
+H.opacity()
+H.transparency()
+
+# # H.animate(0, 'test.mp4')
+
+with h5py.File('h2d.h5', 'w') as f:
+    H.toHdf(f, 'h2d')
+
+with h5py.File('h2d.h5', 'r') as f:
+    H1 = Histogram.fromHdf(f['h2d'])
+
+plt.figure(figsize=(9.5, 5))
+plt.suptitle("Mean, median, and credible interval overlain")
+ax = plt.subplot(121)
+H1.pcolor(cmap='gray_r', colorbar=False)
+H1.plotCredibleIntervals(axis=0)
+H1.plotMedian(axis=0, color='g')
+H1.plotMean(axis=0, color='y')
+plt.legend()
+
+plt.subplot(122, sharex=ax, sharey=ax)
+H1.pcolor(cmap='gray_r', colorbar=False)
+H1.plotCredibleIntervals(axis=1)
+H1.plotMedian(axis=1, color='g')
+H1.plotMean(axis=1, color='y')
+plt.legend()
+
+with h5py.File('h2d.h5', 'w') as f:
+    H.createHdf(f, 'h2d', add_axis=StatArray(np.arange(3.0), name='Easting', units="m"))
+    for i in range(3):
+        H.writeHdf(f, 'h2d', index=i)
+
+with h5py.File('h2d.h5', 'r') as f:
+    H1 = Histogram.fromHdf(f['h2d'], index=0)
+
+plt.figure(figsize=(9.5, 5))
+plt.suptitle("Mean, median, and credible interval overlain")
+ax = plt.subplot(121)
+H1.pcolor(cmap='gray_r', colorbar=False)
+H1.plotCredibleIntervals(axis=0)
+H1.plotMedian(axis=0, color='g')
+H1.plotMean(axis=0, color='y')
+plt.legend()
+
+plt.subplot(122, sharex=ax, sharey=ax)
+H1.pcolor(cmap='gray_r', colorbar=False)
+H1.plotCredibleIntervals(axis=1)
+H1.plotMedian(axis=1, color='g')
+H1.plotMean(axis=1, color='y')
+plt.legend()
+
+with h5py.File('h2d.h5', 'r') as f:
+    H1 = Histogram.fromHdf(f['h2d'])
+
+# H1.pyvista_mesh().save('h3d_read.vtk')
+
+plt.show()

docs/_downloads/09f048ba9debdfb740a9a828a370954b/plot_inference_1d_resolve.py

@@ -0,0 +1,89 @@
+"""
+Running GeoBIPy to invert Resolve data
+++++++++++++++++++++++++++++++++++++++
+"""
+
+import os
+import sys
+import pathlib
+from datetime import timedelta
+import time
+import numpy as np
+from geobipy import Inference3D
+from geobipy import user_parameters
+from geobipy import get_prng
+
+def checkCommandArguments():
+    """Check the users command line arguments. """
+    import argparse
+    # warnings.filterwarnings('error')
+
+    Parser = argparse.ArgumentParser(description="GeoBIPy",
+                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    Parser.add_argument('--index', default=0, type=int, help='job array index 0-18')
+    Parser.add_argument('--data', default=None, help="Data type. Choose from ['skytem_512', 'tempest', 'resolve']")
+    Parser.add_argument('--model', default=None, help="Model type. Choose from ['glacial', 'saline_clay', 'resistive_dolomites', 'resistive_basement', 'coastal_salt_water', 'ice_over_salt_water']")
+
+    return Parser.parse_args()
+
+#%%
+np.random.seed(0)
+
+args = checkCommandArguments()
+sys.path.append(os.getcwd())
+
+models = ['glacial', 'saline_clay', 'resistive_dolomites', 'resistive_basement', 'coastal_salt_water', 'ice_over_salt_water']
+data_type = "Resolve"
+model_type = models[args.index]
+
+#%%
+# The directory where HDF files will be stored
+#%%
+file_path = os.path.join(data_type, model_type)
+pathlib.Path(file_path).mkdir(parents=True, exist_ok=True)
+
+for filename in os.listdir(file_path):
+    try:
+        if os.path.isfile(file_path) or os.path.islink(file_path):
+            os.unlink(file_path)
+    except Exception as e:
+        print('Failed to delete %s. Reason: %s' % (file_path, e))
+
+output_directory = file_path
+
+data_filename = data_type + '_' + model_type
+
+supplementary = "..//..//supplementary//"
+
+parameter_file = supplementary + "//options_files//{}_options".format(data_type)
+inputFile = pathlib.Path(parameter_file)
+assert inputFile.exists(), Exception("Cannot find input file {}".format(inputFile))
+
+output_directory = pathlib.Path(output_directory)
+assert output_directory.exists(), Exception("Make sure the output directory exists {}".format(output_directory))
+
+print('Using user input file {}'.format(parameter_file))
+print('Output files will be produced at {}'.format(output_directory))
+
+kwargs = user_parameters.read(inputFile)
+
+kwargs['n_markov_chains'] = 5000
+
+kwargs['data_filename'] = supplementary + '//data//' + data_filename + '.csv'
+kwargs['system_filename'] = supplementary + "//data//" + kwargs['system_filename']
+
+# Everyone needs the system classes read in early.
+data = kwargs['data_type']._initialize_sequential_reading(kwargs['data_filename'], kwargs['system_filename'])
+
+# Start keeping track of time.
+t0 = time.time()
+
+seed = 146100583096709124601953385843316024947
+prng = get_prng(seed=seed)
+
+inference3d = Inference3D(data, prng=prng)
+inference3d.create_hdf5(directory=output_directory, **kwargs)
+
+print("Created hdf5 files in {} h:m:s".format(str(timedelta(seconds=time.time()-t0))))
+
+inference3d.infer(index=30, **kwargs)