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.
- 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.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.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:
- 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:
- 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.
- burnman.utils.math.l2(x, funca, funcb)[source]¶
Computes the L2 norm for one profile(assumed to be linear between points).
- 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
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.pretty_string_values(popt, pcov, extra_decimal_places=0, combine_value_and_sigma=False)[source]¶
Takes a numpy array of parameters, the corresponding covariance matrix and a set of parameter names and returns the scaled variables and principal 1-s.d.uncertainties (np.sqrt(pcov[i][i])) and scaling factor as three separate lists of strings.
- Parameters:
- Returns:
values, uncertainties and the scaling factors
- Return type:
tuple of 3 lists
- burnman.utils.misc.pretty_print_values(popt, pcov, params, extra_decimal_places=0)[source]¶
Takes a numpy array of parameters, the corresponding covariance matrix and a set of parameter names and prints the scaled variables and principal 1-s.d.uncertainties (np.sqrt(pcov[i][i])) and scaling factor in an easy to read format.
- Parameters:
popt (numpy array) – Parameter values
pcov (2D numpy array) – Variance-covariance matrix
extra_decimal_places (int, optional) – extra precision for values, defaults to 0
- 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.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.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 doneplot_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
orburnman.Layer
.) – Planet or layer to write out to tvel filefilename (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 usingburnman.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:
- 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: