mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
Updated notebooks and added plotDataTypes used for plotting MTdata
This commit is contained in:
GudniRos
@@ -31,26 +31,30 @@
|
||||
},
|
||||
"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",
|
||||
"ct = 20\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,16,1.4)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])\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",
|
||||
"layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)\n",
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"sig_layer1 = .2\n",
|
||||
"sig_layer2 = .2\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
@@ -63,8 +67,7 @@
|
||||
"sigma_0[active] = sig_half\n",
|
||||
"m_0 = np.log(sigma_0[active])\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# Set the mapping# Set the mapping\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"
|
||||
]
|
||||
|
||||
Binary file not shown.
Binary file not shown.
@@ -38,23 +38,23 @@
|
||||
"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",
|
||||
"ct = 20\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,16,1.4)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])\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",
|
||||
"layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)\n",
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"sig_layer1 = .2\n",
|
||||
"sig_layer2 = .2\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
@@ -123,14 +123,14 @@
|
||||
"# 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",
|
||||
"survey.std = np.abs(survey.dobs*std) + 0.01*np.linalg.norm(survey.dobs) #survey.dobs*0 + std\n",
|
||||
"# Assign the data weight\n",
|
||||
"survey.Wd = 1/(abs(survey.dobs)*survey.std)"
|
||||
"survey.Wd = 1/survey.std #(abs(survey.dobs)*survey.std)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -149,13 +149,15 @@
|
||||
"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",
|
||||
"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=mappingExpAct)\n",
|
||||
"reg.alpha_s = 1e-8\n",
|
||||
"reg.smoothModel = True\n",
|
||||
"reg.alpha_s = 1e-9\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"# reg.alpha_xx = 0.001\n",
|
||||
"# Inversion problem\n",
|
||||
"invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)\n",
|
||||
"invProb.counter = C\n",
|
||||
@@ -163,17 +165,17 @@
|
||||
"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",
|
||||
"saveModel.fileName = 'Inversion_NoStoppingregMesh_smoothTrue'\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,saveModel]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -184,46 +186,46 @@
|
||||
"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",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_NoStoppingregMesh_smoothTrue.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",
|
||||
" 0 2.08e+05 2.18e+05 0.00e+00 2.18e+05 4.01e+04 0 \n",
|
||||
" 1 2.08e+05 2.52e+04 2.37e-04 2.52e+04 6.20e+03 0 \n",
|
||||
" 2 2.08e+05 3.46e+03 3.27e-04 3.53e+03 1.03e+03 0 Skip BFGS \n",
|
||||
" 3 2.60e+04 1.72e+03 2.64e-04 1.73e+03 2.61e+02 0 Skip BFGS \n",
|
||||
" 4 2.60e+04 1.03e+03 1.19e-02 1.34e+03 2.44e+02 0 Skip BFGS \n",
|
||||
" 5 2.60e+04 6.85e+02 1.56e-02 1.09e+03 1.35e+02 0 \n",
|
||||
" 6 3.26e+03 5.82e+02 1.78e-02 6.40e+02 1.19e+02 0 \n",
|
||||
" 7 3.26e+03 3.67e+02 4.78e-02 5.22e+02 1.40e+02 0 \n",
|
||||
" 8 3.26e+03 2.60e+02 5.72e-02 4.46e+02 2.06e+02 0 \n",
|
||||
" 9 4.07e+02 2.10e+02 6.04e-02 2.34e+02 6.61e+01 0 \n",
|
||||
" 10 4.07e+02 1.76e+02 9.60e-02 2.15e+02 1.65e+02 0 \n",
|
||||
" 11 4.07e+02 1.56e+02 1.22e-01 2.05e+02 9.68e+01 0 \n",
|
||||
" 12 5.09e+01 1.29e+02 1.28e-01 1.35e+02 5.89e+01 0 \n",
|
||||
" 13 5.09e+01 1.24e+02 1.85e-01 1.34e+02 7.24e+01 0 \n",
|
||||
" 14 5.09e+01 1.19e+02 2.56e-01 1.32e+02 7.95e+01 0 \n",
|
||||
" 15 6.36e+00 1.10e+02 2.10e-01 1.12e+02 5.28e+01 0 \n",
|
||||
" 16 6.36e+00 7.67e+01 3.60e-01 7.89e+01 5.36e+01 0 Skip BFGS \n",
|
||||
" 17 6.36e+00 7.19e+01 3.11e-01 7.39e+01 3.66e+01 0 \n",
|
||||
" 18 7.95e-01 6.99e+01 3.06e-01 7.02e+01 3.89e+01 0 Skip BFGS \n",
|
||||
" 19 7.95e-01 6.92e+01 3.20e-01 6.94e+01 3.72e+01 0 \n",
|
||||
" 20 7.95e-01 6.21e+01 3.14e-01 6.24e+01 8.88e+01 2 Skip BFGS \n",
|
||||
" 21 9.94e-02 5.69e+01 2.69e-01 5.69e+01 3.19e+01 0 \n",
|
||||
" 22 9.94e-02 5.62e+01 2.83e-01 5.63e+01 2.62e+01 0 \n",
|
||||
" 23 9.94e-02 5.46e+01 2.44e-01 5.47e+01 4.84e+01 1 \n",
|
||||
" 24 1.24e-02 5.41e+01 2.08e-01 5.41e+01 4.43e+01 1 Skip BFGS \n",
|
||||
" 25 1.24e-02 5.36e+01 2.09e-01 5.36e+01 4.26e+01 1 \n",
|
||||
" 26 1.24e-02 5.30e+01 1.99e-01 5.30e+01 2.34e+01 0 Skip BFGS \n",
|
||||
" 27 1.55e-03 5.29e+01 1.76e-01 5.29e+01 3.42e+01 0 \n",
|
||||
" 28 1.55e-03 5.28e+01 1.81e-01 5.28e+01 3.56e+01 0 \n",
|
||||
" 29 1.55e-03 5.28e+01 1.75e-01 5.28e+01 3.75e+01 1 \n",
|
||||
" 30 1.94e-04 4.75e+01 1.49e-01 4.75e+01 8.32e+01 1 \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 : |fc-fOld| = 5.2783e+00 <= tolF*(1+|f0|) = 2.1833e+04\n",
|
||||
"0 : |xc-x_last| = 7.4766e+00 <= tolX*(1+|x0|) = 5.1104e+00\n",
|
||||
"0 : |proj(x-g)-x| = 8.3205e+01 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 8.3205e+01 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 30 <= iter = 30\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
@@ -236,14 +238,14 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"modList = []\n",
|
||||
"modFiles = glob('*Inversion_NoStopping.npy')\n",
|
||||
"modFiles = glob('*Inversion_NoStoppingregMesh_smoothTrue.npy')\n",
|
||||
"modFiles.sort()\n",
|
||||
"for f in modFiles:\n",
|
||||
" modList.append(np.load(f))"
|
||||
@@ -251,42 +253,119 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"plt.show()\n"
|
||||
"# simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"# plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# %matplotlib qt\n",
|
||||
"# fig= simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList[-3:-1])\n",
|
||||
"# fig.suptitle('No stopping-useMref ')\n",
|
||||
"# plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"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": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList)\n",
|
||||
"fig = simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])\n",
|
||||
"fig.suptitle('No stopping-useMref')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib qt\n",
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,modList[-3:-1])\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 4.26403069e+01, 3.00486580e+02, 1.66374083e+02,\n",
|
||||
" 7.56749193e+00, 1.95232762e+00, 1.28591697e+01,\n",
|
||||
" 9.25428204e+01, 1.58941963e+02, 1.23328303e+02,\n",
|
||||
" 7.51647545e+01, 4.59620120e+01, 2.84873346e+01,\n",
|
||||
" 1.67854087e+01, 9.65288171e+00, 5.69417961e+00,\n",
|
||||
" 4.11415753e+00, 4.04091270e+00, 5.89702878e+00,\n",
|
||||
" 1.03685013e+01, 1.49615226e+01, 1.78450685e+01,\n",
|
||||
" 2.21298842e+01, 3.14894964e+01, 5.03867813e+01,\n",
|
||||
" 8.57783978e+01, 1.39225191e+02, 2.10683526e+02,\n",
|
||||
" 3.13430053e+02, 4.54348237e+02, 6.35619644e+02,\n",
|
||||
" 8.20077602e+02, 9.86597201e+02, 1.13311026e+03,\n",
|
||||
" 1.22681391e+03, 1.23177188e+03, 1.13462480e+03,\n",
|
||||
" 9.82900015e+02, 7.79370193e+02, 5.53800277e+02,\n",
|
||||
" 3.42690620e+02, 1.91257904e+02, 1.00287868e+02,\n",
|
||||
" 4.40936077e+01, 1.57081219e+01, 4.84819824e+00,\n",
|
||||
" 2.37649055e+00, 2.81483584e+00, 7.91239926e+00,\n",
|
||||
" 3.26193232e+01, 1.10105145e+02, 2.29971251e+02,\n",
|
||||
" 3.48665470e+02, 4.30985817e+02, 4.28578461e+02,\n",
|
||||
" 3.43909988e+02, 2.60243364e+02, 2.15181922e+02,\n",
|
||||
" 2.03527866e+02, 2.38158572e+02, 3.70346010e+02,\n",
|
||||
" 6.61570166e+02, 1.12877904e+03, 1.82064196e+03,\n",
|
||||
" 2.57897162e+03, 3.02323079e+03])"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"1/np.exp(mopt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
|
||||
@@ -53,8 +53,8 @@
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"sig_layer1 = .2\n",
|
||||
"sig_layer2 = .2\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
@@ -102,7 +102,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -130,7 +130,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -154,6 +154,7 @@
|
||||
" reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"# else:\n",
|
||||
"# reg = simpeg.Regularization.Tikhonov(m1d,mapping=mapAct)\n",
|
||||
"reg.smoothModel = False\n",
|
||||
"reg.alpha_s = 1e-8\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"# Inversion problem\n",
|
||||
@@ -163,17 +164,17 @@
|
||||
"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",
|
||||
"saveModel.fileName = 'Inversion_NoStoppingregMesh_smoothFalse'\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,saveModel]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -184,21 +185,48 @@
|
||||
"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",
|
||||
"SimPEG.SaveModelEveryIteration will save your models as: '###-Inversion_NoStoppingregMesh_smoothFalse.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"
|
||||
" 0 7.58e+05 5.84e+06 0.00e+00 5.84e+06 1.45e+06 0 \n",
|
||||
" 1 7.58e+05 6.42e+05 2.33e-04 6.43e+05 1.64e+05 0 \n",
|
||||
" 2 7.58e+05 6.95e+04 7.43e-04 7.00e+04 2.00e+04 0 Skip BFGS \n",
|
||||
" 3 9.47e+04 9.80e+03 1.31e-03 9.92e+03 2.87e+03 0 Skip BFGS \n",
|
||||
" 4 9.47e+04 4.94e+03 5.58e-03 5.47e+03 3.72e+02 0 Skip BFGS \n",
|
||||
" 5 9.47e+04 4.62e+03 7.58e-03 5.33e+03 3.58e+01 0 Skip BFGS \n",
|
||||
" 6 1.18e+04 4.61e+03 7.61e-03 4.70e+03 5.85e+02 0 Skip BFGS \n",
|
||||
" 7 1.18e+04 3.57e+03 3.44e-02 3.97e+03 3.49e+02 0 \n",
|
||||
" 8 1.18e+04 3.27e+03 3.71e-02 3.71e+03 3.65e+02 0 \n",
|
||||
" 9 1.48e+03 3.19e+03 3.67e-02 3.24e+03 5.87e+02 0 \n",
|
||||
" 10 1.48e+03 3.11e+03 5.52e-02 3.19e+03 5.02e+02 0 \n",
|
||||
" 11 1.48e+03 3.10e+03 5.75e-02 3.18e+03 5.14e+02 2 \n",
|
||||
" 12 1.85e+02 3.09e+03 5.75e-02 3.10e+03 5.39e+02 3 \n",
|
||||
" 13 1.85e+02 3.06e+03 7.96e-02 3.07e+03 4.75e+02 0 \n",
|
||||
" 14 1.85e+02 3.06e+03 8.25e-02 3.07e+03 4.80e+02 3 \n",
|
||||
" 15 2.31e+01 3.05e+03 8.50e-02 3.05e+03 4.82e+02 4 \n",
|
||||
" 16 2.31e+01 3.04e+03 1.39e-01 3.04e+03 4.80e+02 2 \n",
|
||||
" 17 2.31e+01 3.01e+03 1.59e-01 3.01e+03 4.72e+02 2 \n",
|
||||
" 18 2.89e+00 3.01e+03 1.26e-01 3.01e+03 4.87e+02 1 \n",
|
||||
" 19 2.89e+00 2.45e+03 1.96e-01 2.45e+03 5.42e+02 1 \n",
|
||||
" 20 2.89e+00 2.11e+03 4.55e-01 2.12e+03 8.59e+02 0 Skip BFGS \n",
|
||||
" 21 3.61e-01 1.95e+03 5.59e-01 1.95e+03 4.88e+02 0 \n",
|
||||
" 22 3.61e-01 1.72e+03 5.09e-01 1.72e+03 3.32e+02 0 \n",
|
||||
" 23 3.61e-01 1.56e+03 6.15e-01 1.56e+03 2.82e+02 0 Skip BFGS \n",
|
||||
" 24 4.52e-02 1.49e+03 7.06e-01 1.49e+03 2.61e+02 1 \n",
|
||||
" 25 4.52e-02 1.45e+03 8.86e-01 1.45e+03 2.72e+02 3 Skip BFGS \n",
|
||||
" 26 4.52e-02 1.41e+03 7.88e-01 1.41e+03 3.08e+02 0 \n",
|
||||
" 27 5.64e-03 1.41e+03 7.02e-01 1.41e+03 2.64e+02 0 \n",
|
||||
" 28 5.64e-03 1.27e+03 7.59e-01 1.27e+03 3.48e+02 0 \n",
|
||||
" 29 5.64e-03 9.93e+02 1.47e+00 9.93e+02 4.30e+02 1 \n",
|
||||
" 30 7.05e-04 8.33e+02 3.08e+00 8.33e+02 5.44e+02 1 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 1.5990e+02 <= tolF*(1+|f0|) = 5.8433e+05\n",
|
||||
"0 : |xc-x_last| = 3.4554e+01 <= tolX*(1+|x0|) = 4.2689e+00\n",
|
||||
"0 : |proj(x-g)-x| = 5.4392e+02 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 5.4392e+02 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 30 <= iter = 30\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -209,7 +237,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -224,7 +252,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -251,13 +279,15 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])\n",
|
||||
"plt.show()\n"
|
||||
"%matplotlib inline\n",
|
||||
"fig = simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])\n",
|
||||
"fig.supplot\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -270,7 +300,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -283,7 +313,29 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"78.694562868192662"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"np.linalg.norm(mopt - reg.mref)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -302,7 +354,7 @@
|
||||
" -6.2146081, -6.2146081, -6.2146081, -6.2146081, -6.2146081])"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -313,37 +365,7 @@
|
||||
},
|
||||
{
|
||||
"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,
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -369,7 +391,7 @@
|
||||
" -5.82076609e-11])"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -414,6 +436,160 @@
|
||||
"m1d.gridN[active]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"self = reg\n",
|
||||
"r1 = self.W * ( self.mapping * (m_0) )\n",
|
||||
"r2 = self.Ws * ( self.mapping * (self.mref) )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"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": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"( self.mapping * (m_0) )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<45x45 sparse matrix of type '<type 'numpy.float64'>'\n",
|
||||
"\twith 45 stored elements in Compressed Sparse Row format>"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"self.Ws"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<SimPEG.Maps.IdentityMap at 0x7f7c7b890390>"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"self.mapping"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"-1.7347234759768071e-18"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"0.5*(r1.dot(r1)-r2.dot(r2))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.029467074824486239"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"r1.dot(r1)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.029467074824486243"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"r2.dot(r2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -65,7 +65,6 @@ def plotMT1DModelData(problem,models,symList=None):
|
||||
|
||||
# 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]
|
||||
|
||||
@@ -0,0 +1,416 @@
|
||||
from matplotlib import pyplot as plt, colors, numpy as np
|
||||
|
||||
|
||||
def rec2nd(structArray):
|
||||
""" Converts a structured/record array to ndarray to do operations on."""
|
||||
return structArray.view((np.float,len(structArray.dtype.names)))
|
||||
|
||||
def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
|
||||
|
||||
indUniFreq = np.where(freq==array['freq'])
|
||||
|
||||
|
||||
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
level = np.logspace(0,-5,31)
|
||||
clevel = np.logspace(0,-4,5)
|
||||
plotNorm = colors.LogNorm()
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
elif colorNorm=='Lin':
|
||||
plotNorm = colors.Normalize()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
return cs
|
||||
|
||||
def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True,mask=None,contourLine=True,useLog=False):
|
||||
|
||||
indUniFreq0 = np.where(freq==arrayList[0]['freq'])
|
||||
indUniFreq1 = np.where(freq==arrayList[1]['freq'])
|
||||
seicmap = plt.get_cmap('RdYlBu')#seismic')
|
||||
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
|
||||
else:
|
||||
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'real':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
|
||||
else:
|
||||
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
|
||||
else:
|
||||
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
cs = ax.tricontourf(x,y,zPlot*100,levels=level*100,cmap=seicmap,extend='both') #,norm=colors.SymLogNorm(1e-2,linscale=2))
|
||||
if contourLine:
|
||||
csl = ax.tricontour(x,y,zPlot*100,levels=clevel*100,colors='k')
|
||||
plt.clabel(csl, fontsize=7, inline=1,fmt='%1.1e',inline_spacing=10)
|
||||
if colorbar:
|
||||
cb = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.1e')
|
||||
for t in cb.ax.get_yticklabels():
|
||||
t.set_rotation(60)
|
||||
t.set_fontsize(8)
|
||||
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
|
||||
def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
|
||||
|
||||
indUniFreq = np.where(freq==array['freq'])
|
||||
|
||||
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
|
||||
if par == 'abs':
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
zPlot = np.abs(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.logspace(-4,0,33)
|
||||
clevel = np.logspace(-4,0,5)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.LogNorm()
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'real':
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
zPlot = np.real(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
|
||||
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'imag':
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
zPlot = np.imag(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
|
||||
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),levels=level,cmap=cmap,norm=plotNorm,edgecolors='k', linewidths=0.5)
|
||||
if colorbar:
|
||||
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
|
||||
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
|
||||
appResFact = 1/(8*np.pi**2*10**(-7))
|
||||
treshold = 1.0 # 1 meter
|
||||
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
|
||||
freq = array['freq'][indUniSta]
|
||||
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indUniSta])
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indUniSta])
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indUniSta])
|
||||
elif par == 'res':
|
||||
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
|
||||
elif par == 'phs':
|
||||
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
|
||||
|
||||
if not pColor:
|
||||
if 'xx' in flag:
|
||||
lab = 'XX'
|
||||
pColor = 'g'
|
||||
elif 'xy' in flag:
|
||||
lab = 'XY'
|
||||
pColor = 'r'
|
||||
elif 'yx' in flag:
|
||||
lab = 'YX'
|
||||
pColor = 'b'
|
||||
elif 'yy' in flag:
|
||||
lab = 'YY'
|
||||
pColor = 'y'
|
||||
|
||||
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
|
||||
|
||||
|
||||
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
|
||||
|
||||
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
x = array['y'][indSect]
|
||||
else:
|
||||
x = array['x'][indSect]
|
||||
y = array['freq'][indSect]
|
||||
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indSect])
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
if cLevel:
|
||||
level = np.logspace(0,-5,31,endpoint=True)
|
||||
clevel = np.logspace(0,-4,5,endpoint=True)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
|
||||
|
||||
elif par == 'ares':
|
||||
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.logspace(zMin,zMax,(zMax-zMin)*8+1,endpoint=True)
|
||||
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
|
||||
plotNorm = colors.LogNorm()
|
||||
|
||||
elif par == 'aphs':
|
||||
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
if cLevel:
|
||||
zMax = cLevel[1]
|
||||
zMin = cLevel[0]
|
||||
else:
|
||||
zMax = (np.ceil(zPlot).max())
|
||||
zMin = (np.floor(zPlot).min())
|
||||
level = np.arange(zMin,zMax+.1,1)
|
||||
clevel = np.arange(zMin,zMax+.1,10)
|
||||
plotNorm = colors.Normalize()
|
||||
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indSect])
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indSect])
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
plotNorm = colors.Normalize()
|
||||
elif colorNorm=='Log':
|
||||
plotNorm = colors.LogNorm()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
# csB.on_mappable_changed(cs)
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
return cs, csB
|
||||
return cs,None
|
||||
|
||||
|
||||
def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=None,contour=True,mask=None,useLog=False):
|
||||
|
||||
def sortInArr(arr):
|
||||
return np.sort(arr,order=['freq','x','y','z'])
|
||||
# Find the index for the slice
|
||||
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
|
||||
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
|
||||
# Extract and sort the mats
|
||||
arr0 = sortInArr(arrayList[0][indSect0])
|
||||
arr1 = sortInArr(arrayList[1][indSect1])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
x0 = arr0['y']
|
||||
x1 = arr1['y']
|
||||
else:
|
||||
x0 = arr0['x']
|
||||
x1 = arr1['x']
|
||||
y0 = arr0['freq']
|
||||
y1 = arr1['freq']
|
||||
|
||||
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
|
||||
else:
|
||||
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
elif par == 'ares':
|
||||
arF = 1/(8*np.pi**2*10**(-7))
|
||||
if useLog:
|
||||
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
|
||||
else:
|
||||
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
|
||||
elif par == 'aphs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
|
||||
else:
|
||||
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
elif par == 'real':
|
||||
if useLog:
|
||||
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
|
||||
else:
|
||||
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('Spectral')
|
||||
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
|
||||
else:
|
||||
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
if cLevel:
|
||||
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
|
||||
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
|
||||
else:
|
||||
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
|
||||
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
|
||||
plotNorm = colors.Normalize()
|
||||
elif colorNorm=='Log':
|
||||
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
|
||||
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
|
||||
plotNorm = colors.LogNorm()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x0,y0,zPlot*100,levels=level*100,cmap=cmap,norm=plotNorm,extend='both')#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x0),np.unique(y0)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.2e')
|
||||
# csB.on_mappable_changed(cs)
|
||||
ax.set_title(flag+' '+par + ' diff',fontsize=8)
|
||||
return cs, csB
|
||||
return cs,None
|
||||
Reference in New Issue
Block a user