Add source code and examples

.gitignore

@@ -0,0 +1,50 @@
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+# MacOS
+.DS_Store
+*/.DS_Store
+
+
+BayesMTGDS
+releases
+tests
+doc
+output/
+*.m
+*.eps
+*.png
+*.jpg
+*.dat
+*.log
+*.logfile

README.md

@@ -1,2 +1,27 @@
-# GEM1DINV-released
-Trans-dimensional Bayesian joint inversion of magnetotelluric and multi-source geomagnetic depth sounding data
+# BayesMTGDS
+
+BayesMTGDS is a standalone C++ code for trans-dimensional Bayesian joint inversion of magnetotelluric (MT) and geomagnetic depth sounding (GDS) responses. Supported data and transfer functions include:
+
+- MT apparent resistivity and phase
+
+- GDS global Qn and Cn responses (a global network of geomagnetic observatory data, satellite data; magnetospheric ring current and ionospheric Sq current (not recommended))
+
+- GDS local Cn responses (single site ground-based data; magnetospheric ring current and ionospheric Sq current (not recommended))
+
+- GDS global-to-local transfer functions (single site ground-based data; magnetospheric ring current and ionospheric Sq current)
+
+### Reference
+The background theory is described in:
+
+- Hongbo Yao, Zhengyong Ren, Jingtian Tang, Rongwen Guo, Jiayong Yan. Trans-dimensional Bayesian joint inversion of magnetotelluric and geomagnetic depth sounding responses to constrain mantle electrical discontinuities. 2022, submitted to Geophysical Journal International
+
+If this helps you, a citation to our paper is appreciated.
+
+### Authors
+
+- Hongbo Yao, Central South University, Email: [email protected]
+
+- Zhengyong Ren, Central South University, Email: [email protected]
+
+### Usage
+Install MPI and type 'make' to compile the source code. Then please see /examples for details.

examples/README.txt

@@ -0,0 +1,13 @@
+/synthetic_data/1D_Kuvshinov
+This example illustrates:
+(1) How to generate synthetic data;
+(2) How to use perform separate and joint inversions.
+
+/synthetic_data/1D_blocky
+This example illustrates how different priors of layer conductivity (uniform distribution, smooth distribution) affect the posterior model distributions
+
+/field_data/TUC
+This folder illustrates a field data example at Tucson geomagnetic observatory
+
+Edit setup.config and type the following command to run the inversion
+> bash run.sh

examples/field_data/TUC/GDS/plot.sh

@@ -0,0 +1,5 @@
+python plot_rms.py
+
+python plot_predicted_Cn_responses.py
+
+python plot_inverse_model.py

examples/field_data/TUC/GDS/plot_inverse_model.py

@@ -0,0 +1,82 @@
+# Copyright (c) 2022 Hongbo Yao
+# Email: [email protected]
+# https://github.com/hongbo-yao
+# https://www.researchgate.net/profile/Hongbo_Yao2 
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['xtick.direction'] = 'in'
+matplotlib.rcParams['ytick.direction'] = 'in'
+
+# plot parameters
+plot_parameters = np.loadtxt('output/plot_parameters.dat')
+z_min = plot_parameters[0]/1000
+z_max = plot_parameters[1]/1000
+
+# mcmc models
+depth_blocks = np.loadtxt('output/depth_blocks.dat')
+depth_blocks = depth_blocks/1000
+sigma_blocks = np.loadtxt('output/sigma_blocks.dat')
+pdf_matrix = np.loadtxt('output/pdf_matrix_normalize.dat')
+[m,n] = pdf_matrix.shape
+pdf_matrix = np.row_stack((pdf_matrix,pdf_matrix[m-1,:]))
+pdf_matrix = np.column_stack((pdf_matrix,pdf_matrix[:,n-1]))
+mean_model = np.loadtxt('output/mean_model.dat')
+mean_model[:,0] = mean_model[:,0]/1000
+median_model = np.loadtxt('output/median_model.dat')
+median_model[:,0] = median_model[:,0]/1000
+credible_interval = np.loadtxt('output/credible_interval.dat')
+credible_interval[:,0] = credible_interval[:,0]/1000
+interface_depth_probability = np.loadtxt('output/interface_depth_probability.dat')
+interface_depth_probability[:,0] = interface_depth_probability[:,0]/1000
+num_of_layer_probability = np.loadtxt('output/num_of_layer_probability.dat')
+
+
+# plot
+linewidth = 1.5
+fontsize = 13
+labelsize = 12
+
+figure = plt.figure(figsize=(13,5))
+
+fig = plt.subplot(1,3,1)
+# cmap: binary, summer
+color = plt.pcolor(sigma_blocks, depth_blocks, pdf_matrix, rasterized=True)
+plt.rcParams['font.size'] = labelsize
+plt.colorbar(color, ax=fig)
+plt.plot(mean_model[:,1], mean_model[:,0], color='m', linewidth=linewidth, label='Mean')
+plt.plot(median_model[:,1], median_model[:,0], color='r', linewidth=linewidth, label='Median')
+plt.plot(credible_interval[:,1], credible_interval[:,0], color='r', linestyle='--', linewidth=linewidth, label='90% credible interval')
+plt.plot(credible_interval[:,2], credible_interval[:,0], color='r', linestyle='--', linewidth=linewidth)
+leg = plt.legend(loc='lower left', facecolor='none', frameon=False, fontsize=fontsize)
+for text in leg.get_texts():
+    text.set_color("white")
+plt.ylim([z_min,z_max])
+fig.invert_yaxis()
+plt.xlabel('Conductivity (log10 S/m)',fontsize=fontsize)
+plt.ylabel('Depth (km)',fontsize=fontsize)
+plt.tick_params(labelsize=labelsize)
+plt.tight_layout()
+
+fig = plt.subplot(1,3,2)
+plt.fill_betweenx(interface_depth_probability[:,0], 0, interface_depth_probability[:,1])
+plt.xlim([0,np.max(interface_depth_probability[:,1])*1.05])
+plt.ylim([z_min,z_max])
+fig.invert_yaxis()
+plt.xlabel('Probability',fontsize=fontsize)
+plt.ylabel('Interface depth (km)',fontsize=fontsize)
+plt.tick_params(labelsize=labelsize)
+plt.tight_layout()
+
+fig = plt.subplot(1,3,3)
+plt.bar(num_of_layer_probability[:,0], num_of_layer_probability[:,1])
+plt.xlabel('Number of layers',fontsize=fontsize)
+plt.ylabel('Probability',fontsize=fontsize)
+plt.tick_params(labelsize=labelsize)
+plt.tight_layout()
+
+plt.savefig('inverse_model.pdf', dpi=600, format='pdf')
+
+plt.show()
+

examples/field_data/TUC/GDS/plot_predicted_Cn_responses.py

@@ -0,0 +1,50 @@
+# Python script for plotting predicted Cn-responses generated by Bayesian1DGEM
+
+# Author:     Hongbo Yao
+# Institute:  School of Geosciences and Info-Physics,
+#             Central South University (CSU)
+# Email:      [email protected]
+# Date:       2022/04/16
+
+# GitHub Page: https://github.com/hongbo-yao
+# Researchgate Page: https://www.researchgate.net/profile/Hongbo_Yao2
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['xtick.direction'] = 'in'
+matplotlib.rcParams['ytick.direction'] = 'in'
+
+# observed data
+observed = np.loadtxt('output/observed_Cn_responses.dat')
+[n_observed,tmp] = observed.shape
+
+# predicted data
+predicted = np.loadtxt('output/predicted_Cn_responses.000000.dat')
+[m,tmp] = predicted.shape
+
+# randomly choosed n_data models for plotting predicted data
+n_data = 50
+loc = np.random.randint(0,m+1,n_data)
+
+# start plotting
+figure = plt.figure(figsize=(6,5))
+fig = plt.subplot(1,1,1)
+plt.errorbar(observed[:,1],observed[:,2],yerr=observed[:,4],fmt='ro',ecolor='r',elinewidth=1,ms=4,mfc='none',mec='r',capthick=1,capsize=3,label='Observed, Re')
+plt.errorbar(observed[:,1],observed[:,3],yerr=observed[:,4],fmt='bo',ecolor='b',elinewidth=1,ms=5,mfc='none',mec='b',capthick=1,capsize=3,label='Observed, Im')
+for i in loc-1:
+   plt.plot(observed[:,1], predicted[i,1:n_observed+1], color='gray', linewidth=0.5)
+   plt.plot(observed[:,1], predicted[i,n_observed+1:2*n_observed+1], color='gray', linewidth=0.5)
+plt.plot(observed[:,1], predicted[loc[n_data-1],1:n_observed+1], color='gray', linewidth=0.5,label='Predicted')
+plt.plot(observed[:,1], predicted[loc[n_data-1],n_observed+1:2*n_observed+1], color='gray', linewidth=0.5)
+plt.xlabel('Period (seconds)',fontsize=14)
+plt.ylabel('C-response (km)',fontsize=14)
+plt.xlim([0.8*observed[0,1],1.2*observed[n_observed-1,1]])
+fig.set_xscale('log') 
+plt.legend(loc='upper left', facecolor='none', frameon=False, fontsize=12)
+plt.tick_params(labelsize=13)
+plt.tight_layout()
+
+plt.savefig('predicted_Cn_responses.pdf', dpi=300, format='pdf')
+
+plt.show()

examples/field_data/TUC/GDS/plot_rms.py

@@ -0,0 +1,39 @@
+# Copyright (c) 2022 Hongbo Yao
+# Email: [email protected]
+# https://github.com/hongbo-yao
+# https://www.researchgate.net/profile/Hongbo_Yao2 
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['xtick.direction'] = 'in'
+matplotlib.rcParams['ytick.direction'] = 'in'
+
+# plot parameters
+plot_parameters = np.loadtxt('output/plot_parameters.dat')
+maxit = plot_parameters[2]
+burnin = plot_parameters[3]
+
+# rms data
+rms = np.loadtxt('output/rms.000000.dat')
+
+# plot
+linewidth = 1.5
+fontsize = 14
+labelsize = 13
+
+figure = plt.figure(figsize=(6,5))
+fig = plt.subplot(1,1,1)
+plt.semilogx(rms[:,0], rms[:,1], '-', linewidth=1)
+# burnin_x = burnin*np.ones(30)
+# burnin_y = np.linspace(0,np.max(rms[:,1]),30)
+# plt.semilogx(burnin_x, burnin_y, 'r--', linewidth=1)
+plt.xlim([1,maxit])
+plt.xlabel('Number of iterations',fontsize=fontsize)
+plt.ylabel('RMS',fontsize=fontsize)
+plt.tick_params(labelsize=labelsize)
+plt.tight_layout()
+
+plt.savefig('RMS.pdf', dpi=600, format='pdf')
+
+plt.show()

examples/field_data/TUC/GDS/run.sh

@@ -0,0 +1,7 @@
+mpirun -np 6 ../../../../BayesMTGDS setup.config | tee BayesMTGDS.log
+
+python plot_rms.py
+
+python plot_predicted_Cn_responses.py
+
+python plot_inverse_model.py

examples/field_data/TUC/GDS/setup.config

@@ -0,0 +1,52 @@
+# -------------------------
+# Observed data:
+# -------------------------
+Observed data file            :  ../field_data/TUC_GDS_data.txt
+
+# -------------------------
+# Common RJMCMC parameters:
+# -------------------------
+Minimum number of layers      :  2        # k_min, default = 2
+Maximum number of layers      :  30       # k_max, default = 40
+Minimum interface depth       :  0        # z_min, default = 0
+Maximum interface depth       :  2000000  # z_max, default = z_cmb
+Minimum log10 conductivity    :  -4       # default = -4
+Maximum log10 conductivity    :  2        # default = 2
+Bitrh probability             :  0.25     # default = 0.25, Bitrh+Death+Move+Update should be 1.0
+Death probability             :  0.25     # default = 0.25
+Move probability              :  0.25     # default = 0.25
+Update probability            :  0.25     # default = 0.25
+Number of mcmc iterations     :  2000000   # default = 500000
+Print every iterations        :  1000     # default = 1000
+
+# -------------------------
+# Priori parameters:
+# -------------------------
+# Minimun layer thicknes        :  1000    # default = (z_max-z_min)/(2*k_max)
+# 0-uniform distribution; 1-normal distribution with reference model constraint; 
+# 2-normal distribution with smooth constraint, default = 2
+Prior distribution type       :  2      
+# Standard derivation for the difference in layer log10 conductivity
+# smaller value will generate smoother models and thiner uncertainty
+Smooth factor                 :  0.5
+
+# -------------------------
+# Proposal parameters:
+# -------------------------
+Birth/death pertubation       :  0.5      # Standard deviation of log10 conductivity for birth/death pertubation, default = (sigma_max-sigma_min)*0.1;
+Update pertubation            :  0.5     # Standard deviation of log10 conductivity for update pertubation, default = (sigma_max-sigma_min)*0.15;
+Move pertubation              :  50000   # Standard deviation of depth in m for move pertubation, default = (z_max-z_min)*0.1;
+
+# -------------------------
+# Post-processing:
+# -------------------------
+Burn in period                :  200000    # default = 10000
+Thin samples                  :  100       # Keep one samples every thin samples, default = 10
+Depth resolution              :  200      # Number of blocks in depth, the larger the higher depth resolution, default = 200
+Conductivity resolution       :  100      # Number of blocks in conductivity, the larger the higher conductivity resolution, default = 100
+Credibl interval              :  0.95     # default = 0.95
+Coarse depth number           :  30       # Number of blocks in depth of coarse model, default = 40, only used when no model_interface_depth.txt
+Write predicted responses     :  1        # 0: do not write responses; 1: only write responses of rank 0; 2: write responses of all processes, default=1
+Write rms misfit              :  2        # 0: do not write rms; 1: only write rms of rank 0; 2: write rms of all processes, default=1
+Output directory              :  output
+Write model samples           :  1

examples/field_data/TUC/Joint/plot.sh

@@ -0,0 +1,7 @@
+python plot_rms.py
+
+python plot_predicted_MT_data.py
+
+python plot_predicted_Cn_responses.py
+
+python plot_inverse_model.py