Skip to content

Commit

Permalink
fix #1 (switching to scientific notation for error graph if param < 0…
Browse files Browse the repository at this point in the history
….01)
  • Loading branch information
sagitta1618 committed Jan 20, 2021
1 parent b054883 commit d26d148
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 22 deletions.
45 changes: 24 additions & 21 deletions src/emagpy/Problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -892,7 +892,7 @@ def solve(obs, pn, spn, alpha, beta, gamma, ini0):
# rrmse = np.sqrt(1/len(app)*np.sum(dataMisfit(cond[:,0], app, ini0)**2)/np.sum(app**2))
# print('ini: RMSE: {:.5f}%'.format(rrmse), ' '.join(
# ['{:.2f}'.format(a) for a in cond[:,0]]))
maxiter = options['maxiter'] if 'maxiter' in options else 1
maxiter = options['maxiter'] if 'maxiter' in options else 3
for l in range(maxiter): # only one iteration as the jacobian doesn't depend on the cond
d = -dataMisfit(cond.flatten(), app, ini0) # NOTE we need minus to get the right direction
LHS = np.dot(J.T, J) + alpha*L
Expand Down Expand Up @@ -2774,6 +2774,8 @@ def resMod2EC(self, fnameECa, fnameResMod, meshType=None, binInt=None, nbins=Non
survey = Survey(fnameECa)
if calib is not None:
survey.gfCorrection(calib=calib)
if meshType is None:
print('Please specify meshType, either "tri" or "quad"')

survey = Survey(fnameECa)
eca = survey.df[['x'] + survey.coils].values
Expand All @@ -2793,9 +2795,10 @@ def resMod2EC(self, fnameECa, fnameResMod, meshType=None, binInt=None, nbins=Non
raise ValueError('You have specified {:d} bins but you only have {:d} data points.'.format(
nbins, eca.shape[0]))

if meshType != 'quad' and 'tri':
print("Please specify mesh type as 'tri' or 'quad'.")
return
if meshType is not None:
if meshType != 'quad' and meshType != 'tri':
print("Please specify mesh type as 'tri' or 'quad'.")
return

#crop ert mesh based on location of EMI calibration data
resmod = resmod[np.where((resmod[:,0] >= minX) & (resmod[:,0] <= maxX)),:][0]
Expand Down Expand Up @@ -2836,24 +2839,24 @@ def resMod2EC(self, fnameECa, fnameResMod, meshType=None, binInt=None, nbins=Non
for i in range(0, len(uniqueX)):
resmod[np.where(resmod[:,0] == uniqueX[i]),1] = newZ #this is done to avoid errors arising from rounding

if meshType == 'tri':
print('Mesh is triangular, topography shift is not implement and will be assumed negligible.')
return
if meshType == 'tri':
print('Mesh is triangular, topography shift is not implement and will be assumed negligible.')
return

#TODO method define upper surface of mesh when topography is involved
x = resmod[:,0]
z = resmod[:,1]
res = resmod[:,2]
xi = np.arange(np.min(x), np.max(x), 0.25)
zi = np.linspace(np.min(z), np.max(z), 15) # 15 layers in Y
xi, zi = np.meshgrid(xi, zi)
resi = griddata((x,z), res, (xi, zi), method='linear') # linear interpolation
x = np.unique(xi)
z = np.unique(zi)
resmodxz = np.array(np.meshgrid(x,z)).T.reshape(-1,2)
res = resi.T.flatten()
resmod = np.concatenate((resmodxz, res[:,None]), axis=1)
resmod = resmod[~np.isnan(resmod[:,2]),:]
#TODO method define upper surface of mesh when topography is involved
x = resmod[:,0]
z = resmod[:,1]
res = resmod[:,2]
xi = np.arange(np.min(x), np.max(x), 0.25)
zi = np.linspace(np.min(z), np.max(z), 15) # 15 layers in Y
xi, zi = np.meshgrid(xi, zi)
resi = griddata((x,z), res, (xi, zi), method='linear') # linear interpolation
x = np.unique(xi)
z = np.unique(zi)
resmodxz = np.array(np.meshgrid(x,z)).T.reshape(-1,2)
res = resi.T.flatten()
resmod = np.concatenate((resmodxz, res[:,None]), axis=1)
resmod = resmod[~np.isnan(resmod[:,2]),:]

midDepths = -np.unique(resmod[:,1])
# compute mean EC for each bin
Expand Down
5 changes: 4 additions & 1 deletion src/emagpy/Survey.py
Original file line number Diff line number Diff line change
Expand Up @@ -860,7 +860,10 @@ def crossOverPointsError(self, coil=None, ax=None, dump=print, minDist=1):
ax.loglog(means, error, '.')
ax.loglog(meansBinned, errorBinned, 'o')
predError = 10**(intercept + slope * np.log10(means))
eq = r'$\epsilon = {:.2f} \times \sigma^{{{:.2f}}}$'.format(10**intercept, slope)
if slope < 0.01 or 10**intercept < 0.01: # switch to scientific notation
eq = r'$\epsilon = {:.2e} \times \sigma^{{{:.2e}}}$'.format(10**intercept, slope)
else:
eq = r'$\epsilon = {:.2f} \times \sigma^{{{:.2f}}}$'.format(10**intercept, slope)
isort = np.argsort(means)
ax.loglog(means[isort], predError[isort], 'k-', label=eq)
ax.legend()
Expand Down

0 comments on commit d26d148

Please sign in to comment.