Skip to content

Commit 200ac5a

Browse files
bobmyhillbobmyhill
committed
sanitise l2 functions
1 parent a5cc064 commit 200ac5a

File tree

11 files changed

+54
-59
lines changed

11 files changed

+54
-59
lines changed

‎burnman/utils/math.py‎

Lines changed: 15 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -346,7 +346,7 @@ def interp_smoothed_array_and_derivatives(
346346
return interps
347347

348348

349-
def compare_l2_norm(depth, calc, obs):
349+
def l2_norm_profiles(depth, calc, obs):
350350
"""
351351
Computes the L2 norm for N profiles at a time (assumed to be linear between points).
352352
@@ -362,12 +362,12 @@ def compare_l2_norm(depth, calc, obs):
362362
"""
363363
err = []
364364
for i in range(len(calc)):
365-
err.append(l2_norm(depth, calc[i], obs[i]))
365+
err.append(l2_norm_profile(depth, calc[i], obs[i]))
366366

367367
return err
368368

369369

370-
def compare_chisqr_factor(calc, obs):
370+
def chisqr_profiles(calc, obs):
371371
"""
372372
Computes the chisqr factor for N profiles at a time. Assumes a 1% a priori uncertainty on the seismic model.
373373
@@ -381,12 +381,12 @@ def compare_chisqr_factor(calc, obs):
381381
"""
382382
err = []
383383
for i in range(len(calc)):
384-
err.append(chisqr_factor(calc[i], obs[i]))
384+
err.append(chisqr_profile(calc[i], obs[i]))
385385

386-
return err
386+
return np.array(err)
387387

388388

389-
def l2_norm(x, funca, funcb):
389+
def l2_norm_profile(x, funca, funcb):
390390
"""
391391
Computes the L2 norm of the difference between two 1D profiles,
392392
assuming linear interpolation between points.
@@ -407,12 +407,11 @@ def l2_norm(x, funca, funcb):
407407
:return: L2 norm (scalar)
408408
:rtype: float
409409
"""
410-
diff = np.array(funca - funcb)
411-
diff = diff * diff
410+
diff = np.square(funca - funcb)
412411
return np.sqrt(integrate.trapezoid(diff, x))
413412

414413

415-
def chisqr_factor(calc, obs):
414+
def chisqr_profile(calc, obs, uncertainties=0.01):
416415
"""
417416
Computes the :math:`\\chi^2` factor for a single profile, assuming a 1 %
418417
uncertainty on the reference (observed) values. This provides a normalized
@@ -428,22 +427,20 @@ def chisqr_factor(calc, obs):
428427
where :math:`N` is the number of elements in the profile.
429428
430429
:param calc: Array of calculated values (e.g., from a model).
431-
:type calc: array-like of float
430+
:type calc: array
432431
433432
:param obs: Array of reference or observed values.
434-
:type obs: array-like of float
433+
:type obs: array
434+
435+
:param uncertainties: Array of uncertainties values.
436+
:type uncertainties: array or float
435437
436438
:return: The :math:`\\chi^2` factor for the profile.
437439
:rtype: float
438440
"""
439441

440-
err = np.empty_like(calc)
441-
for i in range(len(calc)):
442-
err[i] = pow((calc[i] - obs[i]) / (0.01 * np.mean(obs)), 2.0)
443-
444-
err_tot = np.sum(err) / len(err)
445-
446-
return err_tot
442+
err = np.square((calc - obs) / (uncertainties * np.mean(obs)))
443+
return np.sum(err) / len(err)
447444

448445

449446
def independent_row_indices(array):

‎contrib/CHRU2014/paper_onefit.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -325,7 +325,7 @@ def array_to_rock(arr, names):
325325
)
326326

327327
err_vs, err_vphi, err_rho = np.square(
328-
burnman.utils.math.compare_l2_norm(
328+
burnman.utils.math.l2_norm_profiles(
329329
depths / np.mean(depths),
330330
[vs / np.mean(seis_vs), vphi / np.mean(seis_vphi), rho / np.mean(seis_rho)],
331331
[

‎contrib/CHRU2014/paper_opt_pv.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,7 +102,7 @@ def eval_material(amount_perovskite):
102102
def material_error(x):
103103
_, mat_vs, mat_vphi, mat_rho = eval_material(x)
104104
[vs_err, vphi_err, rho_err] = np.square(
105-
burnman.utils.math.compare_l2_norm(
105+
burnman.utils.math.l2_norm_profiles(
106106
depths, [mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
107107
)
108108
)

‎contrib/cider_tutorial_2014/step_2.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ def misfit(phase_1_fraction):
8888
# Here we integrate an L2 difference with depth between our calculated seismic
8989
# profiles and PREM. We then return those misfits.
9090
[vs_err, vphi_err, rho_err] = np.square(
91-
burnman.utils.math.compare_l2_norm(
91+
burnman.utils.math.l2_norm_profiles(
9292
depths, [vs, vphi, density], [seis_vs, seis_vphi, seis_rho]
9393
)
9494
)

‎examples/example_composite_seismic_velocities.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -133,7 +133,7 @@
133133
["density", "v_p", "v_phi", "v_s", "K_S", "G"], seis_p, temperature
134134
)
135135

136-
[vs_err, vphi_err, rho_err] = burnman.utils.math.compare_chisqr_factor(
136+
[vs_err, vphi_err, rho_err] = burnman.utils.math.chisqr_profiles(
137137
[mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
138138
)
139139

‎examples/example_inv_murakami.py‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99

1010
import burnman
1111
from burnman import minerals
12+
from burnman.utils.math import l2_norm_profile
1213

1314
import pymc
1415
import cProfile
@@ -43,12 +44,9 @@ def calc_velocities(amount_pv, iron_pv, iron_fp):
4344
def error(amount_pv, iron_pv, iron_fp):
4445
mat_vp, mat_vs, mat_rho = calc_velocities(amount_pv, iron_pv, iron_fp)
4546

46-
vs_err = np.square(burnman.utils.math.l2_norm(depths, mat_vs, seis_vs)) / 1e9
47-
vp_err = np.square(burnman.utils.math.l2_norm(depths, mat_vp, seis_vp)) / 1e9
48-
den_err = np.square(burnman.utils.math.l2_norm(depths, mat_rho, seis_rho)) / 1e9
49-
50-
# print vs_err, vp_err, den_err
51-
47+
vs_err = l2_norm_profile(depths, mat_vs, seis_vs) / 1.0e3
48+
vp_err = l2_norm_profile(depths, mat_vp, seis_vp) / 1.0e3
49+
den_err = l2_norm_profile(depths, mat_rho, seis_rho) / 1.0e3
5250
return vs_err + vp_err + den_err
5351

5452
# Priors on unknown parameters:
@@ -282,7 +280,7 @@ def theta(amount_pv=amount_pv, iron_pv=iron_pv, iron_fp=iron_fp):
282280
markersize=4,
283281
label="ref",
284282
)
285-
plt.title("density ($\cdot 10^3$ kg/m^$3$)")
283+
plt.title("density ($\\cdot 10^3$ kg/m^$3$)")
286284
plt.legend(loc="upper left")
287285
plt.ylim([4, 8])
288286
plt.savefig("output_figures/example_inv_murakami_2.png")
@@ -525,7 +523,7 @@ def theta(amount_pv=amount_pv, iron_pv=iron_pv, iron_fp=iron_fp):
525523
markersize=4,
526524
label="ref",
527525
)
528-
plt.title("density ($\cdot 10^3$ kg/m$^3$)")
526+
plt.title("density ($\\cdot 10^3$ kg/m$^3$)")
529527
plt.legend(loc="upper left")
530528
plt.ylim([4, 8])
531529
plt.savefig("output_figures/example_inv_murakami_2.png")

‎examples/example_optimize_pv.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
* :class:`burnman.seismic.PREM`
1919
* :func:`burnman.geotherm.brown_shankland`
2020
* :func:`burnman.Material.evaluate`
21-
* :func:`burnman.utils.math.compare_l2_norm`
21+
* :func:`burnman.utils.math.l2_norm_profiles`
2222
2323
*Demonstrates:*
2424
@@ -77,7 +77,7 @@ def material_error(amount_perovskite):
7777
rock.debug_print()
7878
# Calculate errors
7979
[vs_err, vphi_err, rho_err, K_err, G_err] = np.square(
80-
burnman.utils.math.compare_l2_norm(
80+
burnman.utils.math.l2_norm_profiles(
8181
depths,
8282
[mat_vs, mat_vphi, mat_rho, mat_K, mat_G],
8383
[seis_vs, seis_vphi, seis_rho, seis_K, seis_G],

‎examples/example_user_input_material.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -112,7 +112,7 @@ def __init__(self):
112112
["density", "v_s", "v_phi"], seis_p, temperature
113113
)
114114

115-
[vs_err, vphi_err, rho_err] = burnman.utils.math.compare_chisqr_factor(
115+
[vs_err, vphi_err, rho_err] = burnman.utils.math.chisqr_profiles(
116116
[mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
117117
)
118118

‎misc/performance.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77

88
import burnman
99
from burnman import minerals
10-
10+
from burnman.utils.math import l2_norm_profile
1111

1212
if True:
1313

@@ -78,9 +78,9 @@ def calc_velocities(a, b, c):
7878
def error(a, b, c):
7979
mat_vp, mat_vs, mat_rho = calc_velocities(a, b, c)
8080

81-
vs_err = np.square(burnman.utils.math.l2_norm(depths, mat_vs, seis_vs)) / 1e9
82-
vp_err = np.square(burnman.utils.math.l2_norm(depths, mat_vp, seis_vp)) / 1e9
83-
den_err = np.square(burnman.utils.math.l2_norm(depths, mat_rho, seis_rho)) / 1e9
81+
vs_err = l2_norm_profile(depths, mat_vs, seis_vs) / 1.0e3
82+
vp_err = l2_norm_profile(depths, mat_vp, seis_vp) / 1.0e3
83+
den_err = l2_norm_profile(depths, mat_rho, seis_rho) / 1.0e3
8484

8585
return vs_err + vp_err + den_err
8686

‎misc/ref/performance.py.out‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,2 @@
1-
[ 0.03369626 0.00271502 0.06977656 0.00742339 0.00014545]
2-
8130.62642499
1+
[0.03369626 0.00271502 0.06977655 0.00742339 0.00014545]
2+
3183.5272394691788

0 commit comments

Comments
 (0)