mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-18 12:40:30 +08:00
Add notebooks for scipy2015
This commit is contained in:
GudniRos
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because one or more lines are too long
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,321 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the modelling\n",
|
||||
"# Setting up 1D mesh and conductivity models to forward model data.\n",
|
||||
"\n",
|
||||
"# Frequency\n",
|
||||
"nFreq = 31\n",
|
||||
"freqs = np.logspace(3,-3,nFreq)\n",
|
||||
"# Set mesh parameters\n",
|
||||
"ct = 10\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,25,1.3)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,5,-1.2)]),np.ones((3,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],25,-1.3)])\n",
|
||||
"x0 = -np.array([np.sum(np.concatenate((core,bot)))])\n",
|
||||
"# Make the model\n",
|
||||
"m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)\n",
|
||||
"\n",
|
||||
"# Setup model varibles\n",
|
||||
"active = m1d.vectorCCx<0.\n",
|
||||
"layer1 = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-2000.) & (m1d.vectorCCx>=-4000.)\n",
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
"sigma_true[layer1] = sig_layer1\n",
|
||||
"sigma_true[layer2] = sig_layer2\n",
|
||||
"# Extract the model \n",
|
||||
"m_true = np.log(sigma_true[active])\n",
|
||||
"# Make the background model\n",
|
||||
"sigma_0 = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_0[active] = sig_half\n",
|
||||
"m_0 = np.log(sigma_0[active])\n",
|
||||
"\n",
|
||||
"# Set the mapping\n",
|
||||
"actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)\n",
|
||||
"mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the layout of the survey, set the sources and the connected receivers\n",
|
||||
"\n",
|
||||
"# Receivers \n",
|
||||
"rxList = []\n",
|
||||
"for rxType in ['z1dr','z1di']:\n",
|
||||
" rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))\n",
|
||||
"# Source list\n",
|
||||
"srcList =[]\n",
|
||||
"for freq in freqs:\n",
|
||||
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))\n",
|
||||
"# Make the survey\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"survey.mtrue = m_true\n",
|
||||
"# Set the problem\n",
|
||||
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.solver = MumpsSolver\n",
|
||||
"problem.pair(survey)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Read the data\n",
|
||||
"std = 0.05 # 5% std\n",
|
||||
"# Load the files if they exist\n",
|
||||
"if os.path.isfile('MT1D_dtrue.npy') and os.path.isfile('MT1D_dobs.npy'):\n",
|
||||
" d_true = np.load('MT1D_dtrue.npy')\n",
|
||||
" d_obs = np.load('MT1D_dobs.npy')\n",
|
||||
"else:\n",
|
||||
" # Forward model\n",
|
||||
" d_true = survey.dpred(m_true)\n",
|
||||
" np.save('MT1D_dtrue.npy',d_true)\n",
|
||||
" d_obs = d_true + (std*abs(d_true)*np.random.randn(*d_true.shape))\n",
|
||||
" np.save('MT1D_dobs.npy',d_obs)\n",
|
||||
"# Assign the datas to the survey object\n",
|
||||
"survey.dtrue = d_true\n",
|
||||
"survey.dobs = d_obs\n",
|
||||
"survey.std = survey.dobs*0 + std\n",
|
||||
"# Assign the data weight\n",
|
||||
"survey.Wd = 1/(abs(survey.dobs)*survey.std)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the inversion proceedure\n",
|
||||
"\n",
|
||||
"# Define a counter\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"# Set the optimization\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)\n",
|
||||
"opt.counter = C\n",
|
||||
"opt.LSshorten = 0.5\n",
|
||||
"opt.remember('xc')\n",
|
||||
"# Data misfit\n",
|
||||
"dmis = simpeg.DataMisfit.l2_DataMisfit(survey)\n",
|
||||
"# Regularization\n",
|
||||
"# Note: We want you use a mesh the corresponds to the domain we want to solve, the active cells.\n",
|
||||
"if False:\n",
|
||||
" regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)\n",
|
||||
" reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"else:\n",
|
||||
" reg = simpeg.Regularization.Tikhonov(m1d,mapping=mappingExpAct)\n",
|
||||
"reg.alpha_s = 1e-8\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"# Inversion problem\n",
|
||||
"invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)\n",
|
||||
"invProb.counter = C\n",
|
||||
"# Beta cooling\n",
|
||||
"beta = simpeg.Directives.BetaSchedule()\n",
|
||||
"betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"saveModel.fileName = 'Inversion_NoStopping'\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,saveModel]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"SimPEG.InvProblem will set Regularization.mref to m0.\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_NoStopping.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 3.97e+05 1.32e+06 1.77e-07 1.32e+06 3.23e+05 0 \n",
|
||||
" 1 3.97e+05 1.89e+05 3.57e-07 1.89e+05 4.74e+04 0 \n",
|
||||
" 2 3.97e+05 3.88e+04 4.77e-06 3.88e+04 1.11e+04 0 Skip BFGS \n",
|
||||
" 3 4.96e+04 3.69e+04 5.13e-06 3.69e+04 1.05e+04 0 Skip BFGS \n",
|
||||
" 4 4.96e+04 2.71e+04 8.22e-06 2.71e+04 7.80e+03 0 Skip BFGS \n",
|
||||
" 5 4.96e+04 2.34e+04 1.04e-05 2.34e+04 6.80e+03 0 Skip BFGS \n",
|
||||
" 6 6.21e+03 2.11e+04 1.23e-05 2.11e+04 6.19e+03 0 Skip BFGS \n",
|
||||
" 7 6.21e+03 1.27e+04 2.87e-05 1.27e+04 3.95e+03 0 Skip BFGS \n",
|
||||
" 8 6.21e+03 1.07e+04 3.77e-05 1.07e+04 3.41e+03 0 Skip BFGS \n",
|
||||
" 9 7.76e+02 9.53e+03 4.53e-05 9.53e+03 3.09e+03 0 Skip BFGS \n",
|
||||
" 10 7.76e+02 5.51e+03 1.06e-04 5.51e+03 1.97e+03 0 Skip BFGS \n",
|
||||
" 11 7.76e+02 4.57e+03 1.39e-04 4.57e+03 1.69e+03 0 Skip BFGS \n",
|
||||
" 12 9.70e+01 4.02e+03 1.66e-04 4.02e+03 1.53e+03 0 Skip BFGS \n",
|
||||
" 13 9.70e+01 2.27e+03 3.59e-04 2.27e+03 9.73e+02 0 Skip BFGS \n",
|
||||
" 14 9.70e+01 1.83e+03 4.62e-04 1.83e+03 8.24e+02 0 Skip BFGS \n",
|
||||
" 15 1.21e+01 1.57e+03 5.44e-04 1.57e+03 7.36e+02 0 Skip BFGS \n",
|
||||
" 16 1.21e+01 8.91e+02 1.06e-03 8.91e+02 4.57e+02 0 Skip BFGS \n",
|
||||
" 17 1.21e+01 6.95e+02 1.35e-03 6.95e+02 3.66e+02 0 Skip BFGS \n",
|
||||
" 18 1.51e+00 5.92e+02 1.57e-03 5.92e+02 3.13e+02 0 Skip BFGS \n",
|
||||
" 19 1.51e+00 3.56e+02 2.57e-03 3.56e+02 1.75e+02 0 Skip BFGS \n",
|
||||
" 20 1.51e+00 3.22e+02 3.06e-03 3.22e+02 1.53e+02 0 Skip BFGS \n",
|
||||
" 21 1.89e-01 2.75e+02 3.38e-03 2.75e+02 1.25e+02 0 \n",
|
||||
" 22 1.89e-01 1.98e+02 4.78e-03 1.98e+02 8.83e+01 0 Skip BFGS \n",
|
||||
" 23 1.89e-01 1.51e+02 5.79e-03 1.51e+02 7.53e+01 0 \n",
|
||||
" 24 2.37e-02 1.19e+02 6.62e-03 1.19e+02 6.26e+01 0 Skip BFGS \n",
|
||||
" 25 2.37e-02 8.19e+01 1.04e-02 8.19e+01 4.51e+01 0 Skip BFGS \n",
|
||||
" 26 2.37e-02 7.08e+01 1.12e-02 7.08e+01 3.79e+01 1 \n",
|
||||
" 27 2.96e-03 5.49e+01 1.16e-02 5.49e+01 3.51e+01 0 \n",
|
||||
" 28 2.96e-03 4.67e+01 1.23e-02 4.67e+01 3.03e+01 1 \n",
|
||||
" 29 2.96e-03 3.27e+01 1.32e-02 3.27e+01 3.13e+01 0 \n",
|
||||
" 30 3.70e-04 2.51e+01 1.36e-02 2.51e+01 2.30e+01 0 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 7.6530e+00 <= tolF*(1+|f0|) = 1.3164e+05\n",
|
||||
"1 : |xc-x_last| = 4.1960e+00 <= tolX*(1+|x0|) = 4.2689e+00\n",
|
||||
"0 : |proj(x-g)-x| = 2.2988e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 2.2988e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 30 <= iter = 30\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Run the inversion, given the background model as a start.\n",
|
||||
"mopt = inv.run(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"modList = []\n",
|
||||
"modFiles = glob('*Inversion_NoStopping.npy')\n",
|
||||
"modFiles.sort()\n",
|
||||
"for f in modFiles:\n",
|
||||
" modList.append(np.load(f))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList[-3:-1])\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
@@ -0,0 +1,448 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the modelling\n",
|
||||
"# Setting up 1D mesh and conductivity models to forward model data.\n",
|
||||
"\n",
|
||||
"# Frequency\n",
|
||||
"nFreq = 31\n",
|
||||
"freqs = np.logspace(3,-3,nFreq)\n",
|
||||
"# Set mesh parameters\n",
|
||||
"ct = 10\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,25,1.3)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,5,-1.2)]),np.ones((3,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],25,-1.3)])\n",
|
||||
"x0 = -np.array([np.sum(np.concatenate((core,bot)))])\n",
|
||||
"# Make the model\n",
|
||||
"m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)\n",
|
||||
"\n",
|
||||
"# Setup model varibles\n",
|
||||
"active = m1d.vectorCCx<0.\n",
|
||||
"layer1 = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-2000.) & (m1d.vectorCCx>=-4000.)\n",
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
"sigma_true[layer1] = sig_layer1\n",
|
||||
"sigma_true[layer2] = sig_layer2\n",
|
||||
"# Extract the model \n",
|
||||
"m_true = np.log(sigma_true[active])\n",
|
||||
"# Make the background model\n",
|
||||
"sigma_0 = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_0[active] = sig_half\n",
|
||||
"m_0 = np.log(sigma_0[active])\n",
|
||||
"\n",
|
||||
"# Set the mapping\n",
|
||||
"actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)\n",
|
||||
"mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the layout of the survey, set the sources and the connected receivers\n",
|
||||
"\n",
|
||||
"# Receivers \n",
|
||||
"rxList = []\n",
|
||||
"for rxType in ['z1dr','z1di']:\n",
|
||||
" rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))\n",
|
||||
"# Source list\n",
|
||||
"srcList =[]\n",
|
||||
"for freq in freqs:\n",
|
||||
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))\n",
|
||||
"# Make the survey\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"survey.mtrue = m_true\n",
|
||||
"# Set the problem\n",
|
||||
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.solver = MumpsSolver\n",
|
||||
"problem.pair(survey)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Read the data\n",
|
||||
"std = 0.05 # 5% std\n",
|
||||
"# Load the files if they exist\n",
|
||||
"if os.path.isfile('MT1D_dtrue.npy') and os.path.isfile('MT1D_dobs.npy'):\n",
|
||||
" d_true = np.load('MT1D_dtrue.npy')\n",
|
||||
" d_obs = np.load('MT1D_dobs.npy')\n",
|
||||
"else:\n",
|
||||
" # Forward model\n",
|
||||
" d_true = survey.dpred(m_true)\n",
|
||||
" np.save('MT1D_dtrue.npy',d_true)\n",
|
||||
" d_obs = d_true + (std*abs(d_true)*np.random.randn(*d_true.shape))\n",
|
||||
" np.save('MT1D_dobs.npy',d_obs)\n",
|
||||
"# Assign the datas to the survey object\n",
|
||||
"survey.dtrue = d_true\n",
|
||||
"survey.dobs = d_obs\n",
|
||||
"survey.std = np.abs(survey.dobs*std) + 0.01*np.linalg.norm(survey.dobs)\n",
|
||||
"# Assign the data weight\n",
|
||||
"survey.Wd = 1./survey.std"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the inversion proceedure\n",
|
||||
"\n",
|
||||
"# Define a counter\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"# Set the optimization\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)\n",
|
||||
"opt.counter = C\n",
|
||||
"opt.LSshorten = 0.5\n",
|
||||
"opt.remember('xc')\n",
|
||||
"# Data misfit\n",
|
||||
"dmis = simpeg.DataMisfit.l2_DataMisfit(survey)\n",
|
||||
"# Regularization\n",
|
||||
"# Note: We want you use a mesh the corresponds to the domain we want to solve, the active cells.\n",
|
||||
"if True:\n",
|
||||
" regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)\n",
|
||||
" reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"# else:\n",
|
||||
"# reg = simpeg.Regularization.Tikhonov(m1d,mapping=mapAct)\n",
|
||||
"reg.alpha_s = 1e-8\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"# Inversion problem\n",
|
||||
"invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)\n",
|
||||
"invProb.counter = C\n",
|
||||
"# Beta cooling\n",
|
||||
"beta = simpeg.Directives.BetaSchedule()\n",
|
||||
"betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"saveModel.fileName = 'Inversion_NoStoppingregMesh'\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,saveModel]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"SimPEG.InvProblem will set Regularization.mref to m0.\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_NoStoppingregMesh.npy'\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 3.07e+06 5.84e+06 2.95e-02 5.93e+06 1.45e+06 0 \n",
|
||||
" 1 3.07e+06 6.43e+05 2.81e-02 7.29e+05 1.64e+05 0 \n",
|
||||
" 2 3.07e+06 6.97e+04 2.75e-02 1.54e+05 1.99e+04 0 Skip BFGS \n",
|
||||
" 3 3.84e+05 1.01e+04 2.79e-02 2.08e+04 2.95e+03 0 Skip BFGS \n",
|
||||
" 4 3.84e+05 5.52e+03 2.97e-02 1.69e+04 3.93e+02 0 Skip BFGS \n",
|
||||
" 5 3.84e+05 5.44e+03 2.93e-02 1.67e+04 7.12e+01 0 Skip BFGS \n",
|
||||
"------------------------------------------------------------------\n",
|
||||
"0 : ft = 1.6712e+04 <= alp*descent = 1.6712e+04\n",
|
||||
"1 : maxIterLS = 10 <= iterLS = 10\n",
|
||||
"------------------------- End Linesearch -------------------------\n",
|
||||
"The linesearch got broken. Boo.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Run the inversion, given the background model as a start.\n",
|
||||
"mopt = inv.run(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"modList = []\n",
|
||||
"modFiles = glob('*Inversion_NoStopping.npy')\n",
|
||||
"modFiles.sort()\n",
|
||||
"for f in modFiles:\n",
|
||||
" modList.append(np.load(f))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2834: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" if (obj.__array_interface__[\"data\"][0]\n",
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2835: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" != self.__array_interface__[\"data\"][0]):\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])\n",
|
||||
"plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList[-3:-1])\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([-6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081])"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"reg.mref"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([-6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081,\n",
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081])"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m_0"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ -7.62973638e+04, -5.87387576e+04, -4.52321375e+04,\n",
|
||||
" -3.48424297e+04, -2.68503468e+04, -2.07025907e+04,\n",
|
||||
" -1.59735476e+04, -1.23358221e+04, -9.53757167e+03,\n",
|
||||
" -7.38507138e+03, -5.72930192e+03, -4.45563311e+03,\n",
|
||||
" -3.47588787e+03, -2.72223768e+03, -2.14250677e+03,\n",
|
||||
" -1.69655992e+03, -1.35352387e+03, -1.08965000e+03,\n",
|
||||
" -8.86670089e+02, -7.30531699e+02, -6.10425246e+02,\n",
|
||||
" -5.18035666e+02, -4.46966758e+02, -3.92298368e+02,\n",
|
||||
" -3.50245760e+02, -3.17897600e+02, -2.93014400e+02,\n",
|
||||
" -2.68131200e+02, -2.43248000e+02, -2.22512000e+02,\n",
|
||||
" -2.01776000e+02, -1.81040000e+02, -1.63760000e+02,\n",
|
||||
" -1.46480000e+02, -1.29200000e+02, -1.14800000e+02,\n",
|
||||
" -1.00400000e+02, -8.60000000e+01, -7.40000000e+01,\n",
|
||||
" -6.20000000e+01, -5.00000000e+01, -4.00000000e+01,\n",
|
||||
" -3.00000000e+01, -2.00000000e+01, -1.00000000e+01,\n",
|
||||
" -5.82076609e-11])"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"regMesh.gridN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ -7.62973638e+04, -5.87387576e+04, -4.52321375e+04,\n",
|
||||
" -3.48424297e+04, -2.68503468e+04, -2.07025907e+04,\n",
|
||||
" -1.59735476e+04, -1.23358221e+04, -9.53757167e+03,\n",
|
||||
" -7.38507138e+03, -5.72930192e+03, -4.45563311e+03,\n",
|
||||
" -3.47588787e+03, -2.72223768e+03, -2.14250677e+03,\n",
|
||||
" -1.69655992e+03, -1.35352387e+03, -1.08965000e+03,\n",
|
||||
" -8.86670089e+02, -7.30531699e+02, -6.10425246e+02,\n",
|
||||
" -5.18035666e+02, -4.46966758e+02, -3.92298368e+02,\n",
|
||||
" -3.50245760e+02, -3.17897600e+02, -2.93014400e+02,\n",
|
||||
" -2.68131200e+02, -2.43248000e+02, -2.22512000e+02,\n",
|
||||
" -2.01776000e+02, -1.81040000e+02, -1.63760000e+02,\n",
|
||||
" -1.46480000e+02, -1.29200000e+02, -1.14800000e+02,\n",
|
||||
" -1.00400000e+02, -8.60000000e+01, -7.40000000e+01,\n",
|
||||
" -6.20000000e+01, -5.00000000e+01, -4.00000000e+01,\n",
|
||||
" -3.00000000e+01, -2.00000000e+01, -1.00000000e+01])"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m1d.gridN[active]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -0,0 +1,47 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Forward model a data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Define the model\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 182 KiB |
+2
-2
@@ -2,14 +2,14 @@ from simpegEM.FDEM import BaseFDEMProblem
|
||||
from SurveyMT import SurveyMT
|
||||
from DataMT import DataMT
|
||||
from FieldsMT import FieldsMT
|
||||
from SimPEG import SolverLU as SimpegSolver, mkvc
|
||||
from SimPEG import SolverLU as SimpegSolver, Utils, mkvc
|
||||
import numpy as np
|
||||
|
||||
class BaseMTProblem(BaseFDEMProblem):
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
# Set the default pairs of the problem
|
||||
surveyPair = SurveyMT
|
||||
dataPair = DataMT
|
||||
|
||||
@@ -31,7 +31,10 @@ class eForm_psField(BaseMTProblem):
|
||||
@property
|
||||
def sigmaPrimary(self):
|
||||
return self._sigmaPrimary
|
||||
|
||||
@sigmaPrimary.setter
|
||||
def sigmaPrimary(self, val):
|
||||
# Note: TODO add logic for val, make sure it is the correct size.
|
||||
self._sigmaPrimary = val
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
|
||||
@@ -71,7 +71,7 @@ class eForm_ps(BaseMTProblem):
|
||||
|
||||
return C.T*Mmui*C + 1j*omega(freq)*Msig
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
dMe_dsig = self.MeSigmaDeriv( v=u)
|
||||
@@ -94,7 +94,7 @@ class eForm_ps(BaseMTProblem):
|
||||
S_e = Src.S_e(self)
|
||||
return -1j * omega(freq) * S_e
|
||||
|
||||
def getRHSderiv(self, freq, u, v, adjoint=False):
|
||||
def getRHSderiv_m(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
The derivative of the RHS with respect to sigma
|
||||
"""
|
||||
@@ -227,6 +227,7 @@ class eForm_Tp(BaseMTProblem):
|
||||
Abg = C.T*mui*C + 1j*omega(freq)*MeBack
|
||||
|
||||
return Abg*eBG_bp, eBG_bp
|
||||
|
||||
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
|
||||
raise NotImplementedError('getRHSDeriv not implemented yet')
|
||||
return None
|
||||
|
||||
+103
-2
@@ -1,5 +1,7 @@
|
||||
# Utils used for the data,
|
||||
import numpy as np
|
||||
import numpy as np, matplotlib.pyplot as plt, sys
|
||||
import SimPEG as simpeg
|
||||
|
||||
def getAppRes(MTdata):
|
||||
# Make impedance
|
||||
zList = []
|
||||
@@ -21,4 +23,103 @@ def appResPhs(freq,z):
|
||||
return app_res, app_phs
|
||||
|
||||
def rec2ndarr(x,dt=float):
|
||||
return x.view((dt, len(x.dtype.names)))
|
||||
return x.view((dt, len(x.dtype.names)))
|
||||
|
||||
def makeAnalyticSolution(mesh,model,elev,freqs):
|
||||
data1D = []
|
||||
for freq in freqs:
|
||||
anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
|
||||
anaE = anaEd+anaEu
|
||||
anaH = anaHd+anaHu
|
||||
|
||||
anaZ = anaE/anaH
|
||||
# Add to the list
|
||||
data1D.append((freq,0,0,elev,anaZ[0]))
|
||||
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
|
||||
return dataRec
|
||||
|
||||
def plotMT1DModelData(problem,models,symList=None):
|
||||
|
||||
# Setup the figure
|
||||
fontSize = 15
|
||||
|
||||
fig = plt.figure(figsize=[9,7])
|
||||
axM = fig.add_axes([0.075,.1,.25,.875])
|
||||
axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
|
||||
axM.set_xlim(1e-1,1e5)
|
||||
axM.set_ylim(-10000,5000)
|
||||
axM.set_ylabel('Depth [km]',fontsize=fontSize)
|
||||
axR = fig.add_axes([0.42,.575,.5,.4])
|
||||
axR.set_xscale('log')
|
||||
axR.set_yscale('log')
|
||||
axR.invert_xaxis()
|
||||
# axR.set_xlabel('Frequency [Hz]')
|
||||
axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize)
|
||||
|
||||
axP = fig.add_axes([0.42,.1,.5,.4])
|
||||
axP.set_xscale('log')
|
||||
axP.invert_xaxis()
|
||||
axP.set_ylim(0,90)
|
||||
axP.set_xlabel('Frequency [Hz]',fontsize=fontSize)
|
||||
axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize)
|
||||
|
||||
# if not symList:
|
||||
# symList = ['x']*len(models)
|
||||
sys.path.append('/home/gudni/Dropbox/code/python/MTview')
|
||||
import plotDataTypes as pDt
|
||||
# Loop through the models.
|
||||
modelList = [problem.survey.mtrue]
|
||||
modelList.extend(models)
|
||||
if False:
|
||||
modelList = [problem.mapping.sigmaMap*mod for mod in modelList]
|
||||
for nr, model in enumerate(modelList):
|
||||
# Calculate the data
|
||||
if nr==0:
|
||||
data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex')
|
||||
else:
|
||||
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
|
||||
# Plot the data and the model
|
||||
colRat = nr/((len(modelList)-2)*1.)
|
||||
if colRat > 1.:
|
||||
col = 'k'
|
||||
else:
|
||||
col = plt.cm.seismic(1-colRat)
|
||||
# The model - make the pts to plot
|
||||
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
|
||||
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
|
||||
axM.semilogx(modelPts,meshPts,color=col)
|
||||
|
||||
## Data
|
||||
# Appres
|
||||
pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col)
|
||||
# Appphs
|
||||
pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col)
|
||||
try:
|
||||
allData = np.concatenate((allData,simpeg.mkvc(data1D['zyx'],2)),1)
|
||||
except:
|
||||
allData = simpeg.mkvc(data1D['zyx'],2)
|
||||
freq = simpeg.mkvc(data1D['freq'],2)
|
||||
res, phs = appResPhs(freq,allData)
|
||||
|
||||
stdCol = 'gray'
|
||||
axRtw = axR.twinx()
|
||||
axRtw.set_ylabel('Std of log10',color=stdCol)
|
||||
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
|
||||
axPtw = axP.twinx()
|
||||
axPtw.set_ylabel('Std ',color=stdCol)
|
||||
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
|
||||
axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
|
||||
axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
|
||||
|
||||
# Fix labels and ticks
|
||||
|
||||
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
|
||||
axM.set_yticklabels(yMtick)
|
||||
[ l.set_rotation(90) for l in axM.get_yticklabels()]
|
||||
[ l.set_rotation(90) for l in axR.get_yticklabels()]
|
||||
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
|
||||
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
|
||||
for ax in [axM,axR,axP]:
|
||||
ax.xaxis.set_tick_params(labelsize=fontSize)
|
||||
ax.yaxis.set_tick_params(labelsize=fontSize)
|
||||
return fig
|
||||
Reference in New Issue
Block a user