Utilities#

BurnMan has a number of low-level utilities to help achieve common goals. Several of these have already been described in previous sections.

Unit cell#

burnman.utils.unitcell.molar_volume_from_unit_cell_volume(unit_cell_v, z)[source]#

Converts a unit cell volume from Angstroms^3 per unitcell, to m^3/mol.

Parameters:
  • unit_cell_v (float) – Unit cell volumes [A^3/unit cell].

  • z (float) – Number of formula units per unit cell.

Returns:

Volume [m^3/mol]

Return type:

float

burnman.utils.unitcell.cell_parameters_to_vectors(cell_parameters, frame_convention)[source]#

Converts cell parameters to unit cell vectors.

Parameters:
  • cell_parameters (numpy.array (1D)) – An array containing the three lengths of the unit cell vectors [m], and the three angles [degrees]. The first angle (\(\alpha\)) corresponds to the angle between the second and the third cell vectors, the second (\(\beta\)) to the angle between the first and third cell vectors, and the third (\(\gamma\)) to the angle between the first and second vectors.

  • frame_convention (list of three integers) – A list (c) defining the reference frame convention. This function dictates that the c[0]th cell vector is colinear with the c[0]th axis, the c[1]th cell vector is perpendicular to the c[2]th axis, and the c[2]th cell vector is defined to give a right-handed coordinate system. In common crystallographic shorthand, x[c[0]] // a[c[0]], x[c[2]] // a[c[2]]^* (i.e. perpendicular to a[c[0]] and a[c[1]]).

Returns:

The three vectors defining the parallelopiped cell [m].

Return type:

numpy.array (2D)

burnman.utils.unitcell.cell_vectors_to_parameters(vectors, frame_convention)[source]#

Converts unit cell vectors to cell parameters.

Parameters:
  • vectors (numpy.array (2D)) – The three vectors defining the parallelopiped cell [m].

  • frame_convention (list of three integers) – A list (c) defining the reference frame convention. This function dictates that the c[0]th cell vector is colinear with the c[0]th axis, the c[1]th cell vector is perpendicular to the c[2]th axis, and the c[2]th cell vector is defined to give a right-handed coordinate system. In common crystallographic shorthand, x[c[0]] // a[c[0]], x[c[2]] // a[c[2]]^* (i.e. perpendicular to a[c[0]] and a[c[1]]).

Returns:

An array containing the three lengths of the unit cell vectors [m], and the three angles [degrees]. The first angle (\(\alpha\)) corresponds to the angle between the second and the third cell vectors, the second (\(\beta\)) to the angle between the first and third cell vectors, and the third (\(\gamma\)) to the angle between the first and second vectors.

Return type:

numpy.array (1D)

Mathematical#

burnman.utils.math.round_to_n(x, xerr, n)[source]#
burnman.utils.math.unit_normalize(a, order=2, axis=-1)[source]#

Calculates the L2 normalized array of numpy array a of a given order and along a given axis.

burnman.utils.math.float_eq(a, b)[source]#

Test if two floats are almost equal to each other

burnman.utils.math.linear_interpol(x, x1, x2, y1, y2)[source]#

Linearly interpolate to point x, between the points (x1,y1), (x2,y2)

burnman.utils.math.bracket(fn, x0, dx, args=(), ratio=1.618, maxiter=100)[source]#

Given a function and a starting guess, find two inputs for the function that bracket a root.

Parameters:
  • fn (function) – The function to bracket.

  • x0 (float) – The starting guess.

  • dx (float) – Small step for starting the search.

  • args (tuple) – Additional arguments to give to fn.

  • ratio (float) – The step size increases by this ratio every step in the search. Defaults to the golden ratio.

  • maxiter (int) – The maximum number of steps before giving up.

Returns:

xa, xb, fa, fb. xa and xb are the inputs which bracket a root of fn. fa and fb are the values of the function at those points. If the bracket function takes more than maxiter steps, it raises a ValueError.

Return type:

tuple of floats

burnman.utils.math.smooth_array(array, grid_spacing, gaussian_rms_widths, truncate=4.0, mode='inverse_mirror')[source]#

Creates a smoothed array by convolving it with a gaussian filter. Grid resolutions and gaussian RMS widths are required for each of the axes of the numpy array. The smoothing is truncated at a user-defined number of standard deviations. The edges of the array can be padded in a number of different ways given by the ‘mode’ parameter.

Parameters:
  • array (numpy.ndarray) – The array to smooth.

  • grid_spacing (numpy.array of floats) – The spacing of points along each axis.

  • gaussian_rms_widths (numpy.array of floats) – The Gaussian RMS widths/standard deviations for the Gaussian convolution.

  • truncate (float) – The number of standard deviations at which to truncate the smoothing.

  • mode (str) – {‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’, ‘inverse_mirror’} The mode parameter determines how the array borders are handled either by scipy.ndimage.filters.gaussian_filter. Default is ‘inverse_mirror’, which uses burnman.tools.math._pad_ndarray_inverse_mirror().

Returns:

The smoothed array

Return type:

numpy.ndarray

burnman.utils.math.interp_smoothed_array_and_derivatives(array, x_values, y_values, x_stdev=0.0, y_stdev=0.0, truncate=4.0, mode='inverse_mirror', indexing='xy')[source]#

Creates a smoothed array on a regular 2D grid. Smoothing is achieved using burnman.tools.math.smooth_array(). Outputs scipy.interpolate.RegularGridInterpolator() interpolators which can be used to query the array, or its derivatives in the x- and y- directions.

Parameters:
  • array (numpy.array (2D)) – The array to smooth. Each element array[i][j] corresponds to the position x_values[i], y_values[j]

  • x_values (numpy.array (1D)) – The gridded x values over which to create the smoothed grid

  • y_values (numpy.array (1D)) – The gridded y_values over which to create the smoothed grid

  • x_stdev (float) – The standard deviation for the Gaussian filter along the x axis

  • y_stdev (float) – The standard deviation for the Gaussian filter along the y axis

  • truncate (float) – The number of standard deviations at which to truncate the smoothing.

  • mode (str) – {‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’, ‘inverse_mirror’} The mode parameter determines how the array borders are handled either by scipy.ndimage.filters.gaussian_filter. Default is ‘inverse_mirror’, which uses burnman.tools.math._pad_ndarray_inverse_mirror().

  • indexing (str) – {‘xy’, ‘ij’}, optional Cartesian (‘xy’, default) or matrix (‘ij’) indexing of output. See numpy.meshgrid for more details.

Returns:

Three RegularGridInterpolator functors interpolation functions for the smoothed property and the first derivatives with respect to x and y.

Return type:

tuple

burnman.utils.math.compare_l2(depth, calc, obs)[source]#

Computes the L2 norm for N profiles at a time (assumed to be linear between points).

Parameters:
  • depths (array of float) – depths. \([m]\)

  • calc (list of arrays of float) – N arrays calculated values, e.g. [mat_vs,mat_vphi]

  • obs (list of arrays of float) – N arrays of values (observed or calculated) to compare to , e.g. [seis_vs, seis_vphi]

Returns:

array of L2 norms of length N

Return type:

array of floats

burnman.utils.math.compare_chifactor(calc, obs)[source]#

Computes the chi factor for N profiles at a time. Assumes a 1% a priori uncertainty on the seismic model.

Parameters:
  • calc (list of arrays of float) – N arrays calculated values, e.g. [mat_vs,mat_vphi]

  • obs (list of arrays of float) – N arrays of values (observed or calculated) to compare to , e.g. [seis_vs, seis_vphi]

Returns:

error array of length N

Return type:

array of floats

burnman.utils.math.l2(x, funca, funcb)[source]#

Computes the L2 norm for one profile(assumed to be linear between points).

Parameters:
  • x (array of float) – depths \([m]\).

  • funca (list of arrays of float) – array calculated values

  • funcb (list of arrays of float) – array of values (observed or calculated) to compare to

Returns:

L2 norm

Return type:

array of floats

burnman.utils.math.nrmse(x, funca, funcb)[source]#

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

Return type:

array of floats

burnman.utils.math.chi_factor(calc, obs)[source]#

\(\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

Returns:

\(\chi\) factor

Return type:

array of floats

burnman.utils.math.independent_row_indices(array)[source]#

Returns the indices corresponding to an independent set of rows for a given array. The independent rows are determined from the pivots used during row reduction/Gaussian elimination.

Parameters:

array (2D numpy.array of floats) – The input array.

Returns:

The indices corresponding to a set of independent rows of the input array.

Return type:

1D numpy array of integers

burnman.utils.math.array_to_rational_matrix(array)[source]#

Converts a numpy array into a sympy matrix filled with rationals

burnman.utils.math.complete_basis(basis, flip_columns=True)[source]#

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

Parameters:
  • basis (numpy array) – A partial basis

  • flip_columns (bool, optional) – 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

Returns:

basis

Return type:

numpy array

burnman.utils.math.generate_complete_basis(incomplete_basis, array)[source]#

Given a 2D array with independent rows and a second 2D array that spans a larger space, creates a complete basis for the combined array using all the rows of the first array, followed by any required rows of the second array. So, for example, if the first array is: [[1, 0, 0], [1, 1, 0]] and the second array is: [[1, 0, 0], [0, 1, 0], [0, 0, 1]], the complete basis will be: [[1, 0, 0], [1, 1, 0], [0, 0, 1]].

Parameters:
  • incomplete_basis (2D numpy.array) – An array containing the basis to be completed.

  • array (2D numpy.array) – An array spanning the full space for which a basis is required.

Returns:

An array containing the basis vectors spanning both of the input arrays.

Return type:

2D numpy array

burnman.utils.math.is_positive_definite(matrix)[source]#

Checks if a matrix is positive definite

Parameters:

matrix (2D numpy array) – Input matrix

Returns:

Whether or not the matrix is positive definite

Return type:

bool

Miscellaneous#

class burnman.utils.misc.OrderedCounter(iterable=None, /, **kwds)[source]#

Bases: Counter, OrderedDict

Counter that remembers the order elements are first encountered

clear() None.  Remove all items from od.#
copy()#

Return a shallow copy.

elements()#

Iterator over elements repeating each as many times as its count.

>>> c = Counter('ABCABC')
>>> sorted(c.elements())
['A', 'A', 'B', 'B', 'C', 'C']

# Knuth’s example for prime factors of 1836: 2**2 * 3**3 * 17**1 >>> import math >>> prime_factors = Counter({2: 2, 3: 3, 17: 1}) >>> math.prod(prime_factors.elements()) 1836

Note, if an element’s count has been set to zero or is a negative number, elements() will ignore it.

classmethod fromkeys(iterable, v=None)#

Create a new ordered dictionary with keys from iterable and values set to value.

get(key, default=None, /)#

Return the value for key if key is in the dictionary, else default.

items() a set-like object providing a view on D's items#
keys() a set-like object providing a view on D's keys#
most_common(n=None)#

List the n most common elements and their counts from the most common to the least. If n is None, then list all element counts.

>>> Counter('abracadabra').most_common(3)
[('a', 5), ('b', 2), ('r', 2)]
move_to_end(/, key, last=True)#

Move an existing element to the end (or beginning if last is false).

Raise KeyError if the element does not exist.

pop(/, key, default=<unrepresentable>)#

If the key is not found, return the default if given; otherwise, raise a KeyError.

popitem(/, last=True)#

Remove and return a (key, value) pair from the dictionary.

Pairs are returned in LIFO order if last is true or FIFO order if false.

setdefault(/, key, default=None)#

Insert key with a value of default if key is not in the dictionary.

Return the value for key if key is in the dictionary, else default.

subtract(iterable=None, /, **kwds)#

Like dict.update() but subtracts counts instead of replacing them. Counts can be reduced below zero. Both the inputs and outputs are allowed to contain zero and negative counts.

Source can be an iterable, a dictionary, or another Counter instance.

>>> c = Counter('which')
>>> c.subtract('witch')             # subtract elements from another iterable
>>> c.subtract(Counter('watch'))    # subtract elements from another counter
>>> c['h']                          # 2 in which, minus 1 in witch, minus 1 in watch
0
>>> c['w']                          # 1 in which, minus 1 in witch, minus 1 in watch
-1
total()#

Sum of the counts

update(iterable=None, /, **kwds)#

Like dict.update() but add counts instead of replacing them.

Source can be an iterable, a dictionary, or another Counter instance.

>>> c = Counter('which')
>>> c.update('witch')           # add elements from another iterable
>>> d = Counter('watch')
>>> c.update(d)                 # add elements from another counter
>>> c['h']                      # four 'h' in which, witch, and watch
4
values() an object providing a view on D's values#
burnman.utils.misc.copy_documentation(copy_from)[source]#

Decorator @copy_documentation(another_function) will copy the documentation found in a different function (for example from a base class). The docstring applied to some function a() will be

(copied from BaseClass.some_function):
<documentation from BaseClass.some_function>
<optionally the documentation found in a()>
burnman.utils.misc.merge_two_dicts(x, y)[source]#

Given two dicts, merge them into a new dict as a shallow copy.

burnman.utils.misc.flatten(arr)[source]#
burnman.utils.misc.pretty_print_values(popt, pcov, params)[source]#

Takes a numpy array of parameters, the corresponding covariance matrix and a set of parameter names and prints the parameters and principal 1-s.d.uncertainties (np.sqrt(pcov[i][i])) in a nice text based format.

burnman.utils.misc.pretty_print_table(table, use_tabs=False)[source]#

Takes a 2d table and prints it in a nice text based format. If use_tabs=True then only is used as a separator. This is useful for importing the data into other apps (Excel, …). The default is to pad the columns with spaces to make them look neat. The first column is left aligned, while the remainder is right aligned.

burnman.utils.misc.sort_table(table, col=0)[source]#

Sort the table according to the column number

burnman.utils.misc.read_table(filename)[source]#
burnman.utils.misc.array_from_file(filename)[source]#

Generic function to read a file containing floats and commented lines into a 2D numpy array.

Commented lines are prefixed by the characters # or %.

burnman.utils.misc.cut_table(table, min_value, max_value)[source]#
burnman.utils.misc.lookup_and_interpolate(table_x, table_y, x_value)[source]#
burnman.utils.misc.attribute_function(m, attributes, powers=[])[source]#

Function which returns a function which can be used to evaluate material properties at a point. This function allows the user to define the property returned as a string. The function can itself be passed to another function (such as nonlinear_fitting.confidence_prediction_bands()).

Properties can either be simple attributes (e.g. K_T) or a product of attributes, each raised to some power.

Parameters:
  • m (burnman.Material) – The material instance evaluated by the output function.

  • attributes (list of str) – The list of material attributes / properties to be evaluated in the product.

  • powers (list of floats) – The powers to which each attribute should be raised during evaluation.

Returns:

Function which returns the value of product(a_i**p_i) as a function of condition (x = [P, T, V]).

Return type:

function

Tools#

Burnman has a number of high-level tools to help achieve common goals. Several of these have already been described in previous sections.

Plotting#

burnman.tools.plot.pretty_plot()[source]#

Makes pretty plots. Overwrites the matplotlib default settings to allow for better fonts. Slows down plotting

burnman.tools.plot.plot_projected_elastic_properties(mineral, plot_types, axes, n_zenith=31, n_azimuth=91, n_divs=100, n_rots=181, plot_vectors=True)[source]#
Parameters:
  • mineral (burnman.Mineral) – Mineral object on which calculations should be done

  • plot_types (list of str) – Plot types must be one of the following * ‘vp’ - V$_{P}$ (km/s) * ‘vs1’ - ‘V$_{S1}$ (km/s) * ‘vs2’ - V$_{S2}$ (km/s) * ‘vp/vs1’ - V$_{P}$/V$_{S1}$ * ‘vp/vs2’ - V$_{P}$/V$_{S2}$ * ‘s anisotropy’ - S-wave anisotropy (%), 200(vs1s - vs2s)/(vs1s + vs2s)) * ‘linear compressibility’ - Linear compressibility (GPa$^{-1}$) * ‘youngs modulus’ - Youngs Modulus (GPa) * ‘minimum poisson ratio’ - Minimum Poisson ratio * ‘maximum poisson ratio’ - Maximum Poisson ratio

  • axes (matplotlib.pyplot.axes objects) – axes objects to be modified. Must be initialised with projection=’polar’.

  • n_zenith (int) – Number of zeniths (plot resolution).

  • n_azimuth (int) – Number of azimuths (plot resolution).

  • n_divs (int) – Number of divisions for the color scale.

  • n_rots (int) – Number of along-vector rotations to find minimum and maximum Poisson ratios.

  • plot_vectors (bool) – Whether or not to plot vectors for $V_S$ and Poisson axial directions

Output for seismology#

burnman.tools.output_seismo.write_tvel_file(planet_or_layer, modelname='burnmanmodel', background_model=None)[source]#

Function to write input file for obspy travel time calculations. Note: Because density isn’t defined for most 1D seismic models, densities are output as zeroes. The tvel format has a column for density, but this column is not used by obspy for travel time calculations.

Parameters:
  • planet_or_layer (burnman.Planet or burnman.Layer.) – Planet or layer to write out to tvel file

  • filename (str) – Filename to read to.

  • background_model (burnman.seismic.Seismic1DModel) – 1D seismic model to fill in parts of planet (likely to be an earth model) that aren’t defined by layer (only need when using burnman.Layer)

burnman.tools.output_seismo.write_axisem_input(layers, modelname='burnmanmodel_foraxisem', axisem_ref='axisem_prem_ani_noocean.txt', plotting=False)[source]#

Writing velocities and densities to AXISEM (www.axisem.info) input file. The input can be a single layer, or a list of layers taken from a planet (planet.layers). Currently this function will implement explicit discontinuities between layers in the seismic model. Currently this function is only set for Earth.

Parameters:
  • layers (list of one or more burnman.Layer) – List of layers to put in AXISEM file.

  • modelname (str) – Name of model, appears in name of output file.

  • axisem_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

  • plotting (bool) – Choose whether to show plot of the old model and replaced model.

burnman.tools.output_seismo.write_mineos_input(layers, modelname='burnmanmodel_formineos', mineos_ref='mineos_prem_noocean.txt', plotting=False)[source]#

Writing velocities and densities to Mineos (https://geodynamics.org/cig/software/mineos/) input file Note: currently, this function only honors the discontinuities already in the synthetic input file, so it is best to only replace certain layers with burnman values

Parameters:
  • layers (list of one or more burnman.Layer) – List of layers to put in Mineos file.

  • modelname (str) – Name of model, appears in name of output file.

  • mineos_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

  • plotting (bool) – Choose whether to show plot of the old model and replaced model.

Equations of state#

burnman.tools.eos.check_eos_consistency(m, P=1000000000.0, T=300.0, tol=0.0001, verbose=False, including_shear_properties=True, equilibration_function=None)[source]#

Checks that numerical derivatives of the Gibbs energy of a mineral under given conditions are equal to those provided analytically by the equation of state.

Parameters:
  • m (burnman.Mineral) – The mineral for which the equation of state is to be checked for consistency.

  • P (float) – The pressure at which to check consistency.

  • T (float) – The temperature at which to check consistency.

  • tol (float) – The fractional tolerance for each of the checks.

  • verbose (bool) – Decide whether to print information about each check.

  • including_shear_properties (bool) – Decide whether to check shear information, which is pointless for liquids and equations of state without shear modulus parameterizations.

  • equilibration_function (function) – Function to internally equilibrate object. Called after every set_state. Takes the mineral object as its only argument.

Returns:

Boolean stating whether all checks have passed.

Return type:

bool

burnman.tools.eos.check_anisotropic_eos_consistency(m, P=1000000000.0, T=300.0, tol=0.0001, verbose=False, equilibration_function=None)[source]#

Checks that numerical derivatives of the Gibbs energy of an anisotropic mineral under given conditions are equal to those provided analytically by the equation of state.

Parameters:
  • m (burnman.AnisotropicMineral) – The anisotropic mineral for which the equation of state is to be checked for consistency.

  • P (float) – The pressure at which to check consistency.

  • T (float) – The temperature at which to check consistency.

  • tol (float) – The fractional tolerance for each of the checks.

  • verbose (bool) – Decide whether to print information about each check.

  • equilibration_function (function) – Function to internally equilibrate object. Called after every set_state. Takes the mineral object as its only argument.

Returns:

Boolean stating whether all checks have passed.

Return type:

bool