Merge pull request #653 from bobmyhill/fix_l2_plus_tests

Fix l2 plus tests

Author: bobmyhill

burnman/utils/math.py

@@ -30,13 +30,6 @@ def unit_normalize(a, order=2, axis=-1):
     return a / np.expand_dims(l2, axis)[0][0]
 
 
-def float_eq(a, b):
-    """
-    Test if two floats are almost equal to each other
-    """
-    return abs(a - b) < 1e-10 * max(1e-5, abs(a), abs(b))
-
-
 def linear_interpol(x, x1, x2, y1, y2):
     """
     Linearly interpolate to point x, between
@@ -353,7 +346,7 @@ def interp_smoothed_array_and_derivatives(
     return interps
 
 
-def compare_l2(depth, calc, obs):
+def l2_norm_profiles(depth, calc, obs):
     """
     Computes the L2 norm for N profiles at a time (assumed to be linear between points).
 
@@ -369,91 +362,85 @@ def compare_l2(depth, calc, obs):
     """
     err = []
     for i in range(len(calc)):
-        err.append(l2(depth, calc[i], obs[i]))
+        err.append(l2_norm_profile(depth, calc[i], obs[i]))
 
     return err
 
 
-def compare_chifactor(calc, obs):
+def chisqr_profiles(calc, obs):
     """
-    Computes the chi factor for N profiles at a time. Assumes a 1% a priori uncertainty on the seismic model.
-
+    Computes the chisqr factor for N profiles at a time. Assumes a 1% a priori uncertainty on the seismic model.
 
     :type calc: list of arrays of float
     :param calc: N arrays calculated values, e.g. [mat_vs,mat_vphi]
     :type obs: list of arrays of float
-    :param obs: N arrays of values (observed or calculated) to compare to , e.g. [seis_vs, seis_vphi]
+    :param obs: N arrays of values (observed or calculated) to compare to, e.g. [seis_vs, seis_vphi]
 
     :returns: error array of length N
     :rtype: array of floats
     """
     err = []
     for i in range(len(calc)):
-        err.append(chi_factor(calc[i], obs[i]))
+        err.append(chisqr_profile(calc[i], obs[i]))
 
-    return err
+    return np.array(err)
 
 
-def l2(x, funca, funcb):
+def l2_norm_profile(x, funca, funcb):
     """
-    Computes the L2 norm for one profile(assumed to be linear between points).
+    Computes the L2 norm of the difference between two 1D profiles,
+    assuming linear interpolation between points.
 
-    :type x: array of float
-    :param x: depths :math:`[m]`.
-    :type funca: list of arrays of float
-    :param funca: array calculated values
-    :type funcb: list of arrays of float
-    :param funcb: array of values (observed or calculated) to compare to
+    This is equivalent to the square root of the integrated squared
+    difference between the two functions over the given depth range.
+    The integration is performed using the trapezoidal rule.
 
-    :returns: L2 norm
-    :rtype: array of floats
-    """
-    diff = np.array(funca - funcb)
-    diff = diff * diff
+    :param x: Array of depth values [m], must be 1D and monotonically increasing.
+    :type x: array-like of float
 
-    return integrate.trapezoid(diff, x)
+    :param funca: First profile (e.g., model or predicted values), must be same length as x.
+    :type funca: array-like of float
 
+    :param funcb: Second profile (e.g., observed or reference values), same shape as funca.
+    :type funcb: array-like of float
 
-def nrmse(x, funca, funcb):
+    :return: L2 norm (scalar)
+    :rtype: float
     """
-    Normalized root mean square error for one profile
-    :type x: array of float
-    :param x: depths in m.
-    :type funca: list of arrays of float
-    :param funca: array calculated values
-    :type funcb: list of arrays of float
-    :param funcb: array of values (observed or calculated) to compare to
-
-    :returns: RMS error
-    :rtype: array of floats
+    diff = np.square(funca - funcb)
+    return np.sqrt(integrate.trapezoid(diff, x))
+
+
+def chisqr_profile(calc, obs, uncertainties=0.01):
     """
-    diff = np.array(funca - funcb)
-    diff = diff * diff
-    rmse = np.sqrt(np.sum(diff) / x)
-    nrmse = rmse / (np.max(funca) - np.min(funca))
+    Computes the :math:`\\chi^2` factor for a single profile, assuming a 1 %
+    uncertainty on the reference (observed) values. This provides a normalized
+    measure of the deviation between calculated and reference values.
 
-    return nrmse
+    The :math:`\\chi^2` factor is calculated as:
 
+    .. math::
 
-def chi_factor(calc, obs):
-    """
-    :math:`\\chi` factor for one profile assuming 1% uncertainty on the reference model (obs)
-    :type calc: list of arrays of float
-    :param calc: array calculated values
-    :type obs: list of arrays of float
-    :param obs: array of reference values to compare to
+        \\chi^2 = \\frac{1}{N} \\sum_{i=1}^{N}
+                  \\left( \\frac{\\text{calc}_i - \\text{obs}_i}{0.01 \\cdot \\text{mean}(\\text{obs})} \\right)^2
 
-    :returns: :math:`\\chi` factor
-    :rtype: array of floats
-    """
+    where :math:`N` is the number of elements in the profile.
 
-    err = np.empty_like(calc)
-    for i in range(len(calc)):
-        err[i] = pow((calc[i] - obs[i]) / (0.01 * np.mean(obs)), 2.0)
+    :param calc: Array of calculated values (e.g., from a model).
+    :type calc: array
 
-    err_tot = np.sum(err) / len(err)
+    :param obs: Array of reference or observed values.
+    :type obs: array
+
+    :param uncertainties: Array of uncertainties values.
+    :type uncertainties: array or float
+
+    :return: The :math:`\\chi^2` factor for the profile.
+    :rtype: float
+    """
 
-    return err_tot
+    err = np.square((calc - obs) / (uncertainties * np.mean(obs)))
+    return np.sum(err) / len(err)
 
 
 def independent_row_indices(array):
@@ -495,27 +482,33 @@ def array_to_rational_matrix(array):
 
 def complete_basis(basis, flip_columns=True):
     """
-    Creates a full basis by filling remaining rows with
-    selected rows of the identity matrix that have indices
-    not in the column pivot list of the basis RREF
-
-    :param basis: A partial basis
-    :type basis: numpy array
-    :param flip_columns: whether to calculate the RREF
-        from a flipped basis. This can be useful, for example,
-        when dependent columns are known to be further to the
-        right in the basis. defaults to True
+    Extends a partial basis (with fewer rows than columns) to a full basis
+    by adding rows from the identity matrix that are linearly independent
+    of the existing rows.
+
+    This is done by computing the reduced row echelon form (RREF) of the basis
+    to identify pivot columns. Identity rows corresponding to non-pivot columns
+    are added to complete the basis.
+
+    :param basis: A 2D NumPy array representing a partial basis, with shape (n, m)
+                  where n < m.
+    :type basis: numpy.ndarray
+
+    :param flip_columns: If True, the basis columns are flipped (reversed) before
+                         computing the RREF. This is helpful if the dependent columns
+                         are expected to appear later in the matrix. Default is True.
     :type flip_columns: bool, optional
-    :return: basis
-    :rtype: numpy array
-    """
 
+    :return: A completed basis as a NumPy array with shape (m, m), containing the
+             original rows plus additional rows from the identity matrix.
+    :rtype: numpy.ndarray
+    """
     n, m = basis.shape
     if n < m:
         if flip_columns:
             pivots = list(m - 1 - np.array(Matrix(np.flip(basis, axis=1)).rref()[1]))
         else:
-            pivots = list(Matrix(basis)).rref()[1]
+            pivots = list(Matrix(basis).rref())[1]
 
         return np.concatenate(
             (basis, np.identity(m)[[i for i in range(m) if i not in pivots], :]), axis=0

contrib/CHRU2014/paper_onefit.py

@@ -324,14 +324,16 @@ def array_to_rock(arr, names):
         ["rho", "v_p", "v_s", "v_phi", "K_S", "G"], pressure, temperature
     )
 
-    err_vs, err_vphi, err_rho = burnman.utils.math.compare_l2(
-        depths / np.mean(depths),
-        [vs / np.mean(seis_vs), vphi / np.mean(seis_vphi), rho / np.mean(seis_rho)],
-        [
-            seis_vs / np.mean(seis_vs),
-            seis_vphi / np.mean(seis_vphi),
-            seis_rho / np.mean(seis_rho),
-        ],
+    err_vs, err_vphi, err_rho = np.square(
+        burnman.utils.math.l2_norm_profiles(
+            depths / np.mean(depths),
+            [vs / np.mean(seis_vs), vphi / np.mean(seis_vphi), rho / np.mean(seis_rho)],
+            [
+                seis_vs / np.mean(seis_vs),
+                seis_vphi / np.mean(seis_vphi),
+                seis_rho / np.mean(seis_rho),
+            ],
+        )
     )
     error = np.sum([err_rho, err_vphi, err_vs])
 

contrib/CHRU2014/paper_opt_pv.py

@@ -97,14 +97,14 @@ def eval_material(amount_perovskite):
         mat_rho, mat_vs, mat_vphi = rock.evaluate(
             ["rho", "v_s", "v_phi"], seis_p, temperature
         )
-        # [rho_err,vphi_err,vs_err]=burnman.utils.compare_chifactor(mat_vs,mat_vphi,mat_rho,seis_vs,seis_vphi,seis_rho)
-
         return seis_p, mat_vs, mat_vphi, mat_rho
 
     def material_error(x):
         _, mat_vs, mat_vphi, mat_rho = eval_material(x)
-        [vs_err, vphi_err, rho_err] = burnman.utils.math.compare_l2(
-            depths, [mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
+        [vs_err, vphi_err, rho_err] = np.square(
+            burnman.utils.math.l2_norm_profiles(
+                depths, [mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
+            )
         )
         scale = 2700e3 - 850e3
         return vs_err / scale, vphi_err / scale

contrib/cider_tutorial_2014/step_2.py

@@ -87,8 +87,10 @@ def misfit(phase_1_fraction):
 
         # Here we integrate an L2 difference with depth between our calculated seismic
         # profiles and PREM.  We then return those misfits.
-        [vs_err, vphi_err, rho_err] = burnman.utils.math.compare_l2(
-            depths, [vs, vphi, density], [seis_vs, seis_vphi, seis_rho]
+        [vs_err, vphi_err, rho_err] = np.square(
+            burnman.utils.math.l2_norm_profiles(
+                depths, [vs, vphi, density], [seis_vs, seis_vphi, seis_rho]
+            )
         )
 
         return vs_err, vphi_err, rho_err

examples/example_composite_seismic_velocities.py

@@ -133,7 +133,7 @@
         ["density", "v_p", "v_phi", "v_s", "K_S", "G"], seis_p, temperature
     )
 
-    [vs_err, vphi_err, rho_err] = burnman.utils.math.compare_chifactor(
+    [vs_err, vphi_err, rho_err] = burnman.utils.math.chisqr_profiles(
         [mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
     )
 

examples/example_inv_murakami.py

@@ -9,6 +9,7 @@
 
 import burnman
 from burnman import minerals
+from burnman.utils.math import l2_norm_profile
 
 import pymc
 import cProfile
@@ -43,12 +44,9 @@ def calc_velocities(amount_pv, iron_pv, iron_fp):
     def error(amount_pv, iron_pv, iron_fp):
         mat_vp, mat_vs, mat_rho = calc_velocities(amount_pv, iron_pv, iron_fp)
 
-        vs_err = burnman.utils.math.l2(depths, mat_vs, seis_vs) / 1e9
-        vp_err = burnman.utils.math.l2(depths, mat_vp, seis_vp) / 1e9
-        den_err = burnman.utils.math.l2(depths, mat_rho, seis_rho) / 1e9
-
-        # print vs_err, vp_err, den_err
-
+        vs_err = l2_norm_profile(depths, mat_vs, seis_vs) / 1.0e3
+        vp_err = l2_norm_profile(depths, mat_vp, seis_vp) / 1.0e3
+        den_err = l2_norm_profile(depths, mat_rho, seis_rho) / 1.0e3
         return vs_err + vp_err + den_err
 
     # Priors on unknown parameters:
@@ -282,7 +280,7 @@ def theta(amount_pv=amount_pv, iron_pv=iron_pv, iron_fp=iron_fp):
             markersize=4,
             label="ref",
         )
-        plt.title("density ($\cdot 10^3$ kg/m^$3$)")
+        plt.title("density ($\\cdot 10^3$ kg/m^$3$)")
         plt.legend(loc="upper left")
         plt.ylim([4, 8])
         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):
             markersize=4,
             label="ref",
         )
-        plt.title("density ($\cdot 10^3$ kg/m$^3$)")
+        plt.title("density ($\\cdot 10^3$ kg/m$^3$)")
         plt.legend(loc="upper left")
         plt.ylim([4, 8])
         plt.savefig("output_figures/example_inv_murakami_2.png")

examples/example_optimize_pv.py

@@ -18,7 +18,7 @@
 * :class:`burnman.seismic.PREM`
 * :func:`burnman.geotherm.brown_shankland`
 * :func:`burnman.Material.evaluate`
-* :func:`burnman.utils.math.compare_l2`
+* :func:`burnman.utils.math.l2_norm_profiles`
 
 *Demonstrates:*
 
@@ -76,10 +76,12 @@ def material_error(amount_perovskite):
         print("Calculations are done for:")
         rock.debug_print()
         # Calculate errors
-        [vs_err, vphi_err, rho_err, K_err, G_err] = burnman.utils.math.compare_l2(
-            depths,
-            [mat_vs, mat_vphi, mat_rho, mat_K, mat_G],
-            [seis_vs, seis_vphi, seis_rho, seis_K, seis_G],
+        [vs_err, vphi_err, rho_err, K_err, G_err] = np.square(
+            burnman.utils.math.l2_norm_profiles(
+                depths,
+                [mat_vs, mat_vphi, mat_rho, mat_K, mat_G],
+                [seis_vs, seis_vphi, seis_rho, seis_K, seis_G],
+            )
         )
         # Normalize errors
         vs_err = vs_err / np.mean(seis_vs) ** 2.0

examples/example_user_input_material.py

@@ -112,7 +112,7 @@ def __init__(self):
         ["density", "v_s", "v_phi"], seis_p, temperature
     )
 
-    [vs_err, vphi_err, rho_err] = burnman.utils.math.compare_chifactor(
+    [vs_err, vphi_err, rho_err] = burnman.utils.math.chisqr_profiles(
         [mat_vs, mat_vphi, mat_rho], [seis_vs, seis_vphi, seis_rho]
     )
 

misc/performance.py

@@ -1,12 +1,13 @@
-# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
-# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
+# for the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
 import numpy as np
 
 import burnman
 from burnman import minerals
-
+from burnman.utils.math import l2_norm_profile
 
 if True:
 
@@ -77,9 +78,9 @@ def calc_velocities(a, b, c):
     def error(a, b, c):
         mat_vp, mat_vs, mat_rho = calc_velocities(a, b, c)
 
-        vs_err = burnman.utils.math.l2(depths, mat_vs, seis_vs) / 1e9
-        vp_err = burnman.utils.math.l2(depths, mat_vp, seis_vp) / 1e9
-        den_err = burnman.utils.math.l2(depths, mat_rho, seis_rho) / 1e9
+        vs_err = l2_norm_profile(depths, mat_vs, seis_vs) / 1.0e3
+        vp_err = l2_norm_profile(depths, mat_vp, seis_vp) / 1.0e3
+        den_err = l2_norm_profile(depths, mat_rho, seis_rho) / 1.0e3
 
         return vs_err + vp_err + den_err
 

misc/ref/performance.py.out

@@ -1,2 +1,2 @@
-[ 0.03369626  0.00271502  0.06977656  0.00742339  0.00014545]
-8130.62642499
+[0.03369626 0.00271502 0.06977655 0.00742339 0.00014545]
+3183.5272394691788