Optimization

class burnman.optimize.composition_fitting.DummyCompositionSolution(endmember_element_formulae, endmember_site_formulae)

Bases: Solution

This is a dummy base class for a solution object to facilitate composition fitting when no solution model has been prepared. The model is initialized with appropriate chemical formulae for each endmember, and can do all basic compositional processing that doesn’t involve any material properties.

Parameters:

• endmember_element_formulae (list of str) – Formulae for each of the independent endmembers. e.g. [‘Mg2SiO4’, ‘Fe2SiO4’].

• endmember_site_formulae (list of str) – Site formulae for each of the independent endmembers, in the same order as endmember_element_formulae. e.g. [‘[Mg]2SiO4’, ‘[Fe]2SiO4’].

property endmember_formulae

A list of formulae for all the endmember in the solution.

property C_p

Alias for molar_heat_capacity_p()

property C_v

Alias for molar_heat_capacity_v()

property G

Alias for shear_modulus()

property H

Alias for molar_enthalpy()

property K_S

Alias for isentropic_bulk_modulus_reuss()

property K_T

Alias for isothermal_bulk_modulus_reuss()

property P

Alias for pressure()

property S

Alias for molar_entropy()

property T

Alias for temperature()

property V

Alias for molar_volume()

property activities

Returns a list of endmember activities [unitless].

property activity_coefficients

Returns a list of endmember activity coefficients (gamma = activity / ideal activity) [unitless].

property alpha

Alias for thermal_expansivity()

property beta_S

Alias for isentropic_compressibility_reuss()

property beta_T

Alias for isothermal_compressibility_reuss()

property bulk_sound_velocity

Returns bulk sound speed of the solution [m/s]. Aliased with self.v_phi.

property compositional_basis

_summary_

Returns:

_description_

Return type:

_type_

property compositional_null_basis

An array N such that N.b = 0 for all bulk compositions that can be produced with a linear sum of the endmembers in the solution.

copy()

debug_print(indent='')

Print a human-readable representation of this Material.

property density

Returns density of the solution [kg/m^3]. Aliased with self.rho.

property dependent_element_indices

The element indices not included in the independent list.

property elements

A list of the elements which could be contained in the solution, returned in the IUPAC element order.

property endmember_names

A list of names for all the endmember in the solution.

property endmembers

property energy

Alias for molar_internal_energy()

property entropy_hessian

Returns an array containing the second compositional derivative of the entropy [J/K]. Property specific to solutions.

evaluate(vars_list, pressures, temperatures, molar_fractions=None)

Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.

Parameters:

• vars_list (list of strings) – Variables to be returned for given conditions

• pressures (numpy.array, n-dimensional) – ndlist or ndarray of float of pressures in [Pa].

• temperatures (numpy.array, n-dimensional) – ndlist or ndarray of float of temperatures in [K].

Returns:

List or array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] for temperatures[i] and pressures[i]. Attempts to return an array, falls back to a list if the returned properties have different shapes.

Return type:

list or numpy.array, n-dimensional

evaluate_with_volumes(vars_list, volumes, temperatures, molar_fractions=None)

Returns an array of material properties requested through a list of strings at given volume and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.

Parameters:

• vars_list (list of strings) – Variables to be returned for given conditions

• volumes (numpy.array, n-dimensional) – ndlist or ndarray of float of volumes in [m^3].

• temperatures (numpy.array, n-dimensional) – ndlist or ndarray of float of temperatures in [K].

Returns:

List or array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] for temperatures[i] and pressures[i]. Attempts to return an array, falls back to a list if the returned properties have different shapes.

Return type:

list or numpy.array, n-dimensional

property excess_enthalpy

Returns excess molar enthalpy [J/mol]. Property specific to solutions.

property excess_entropy

Returns excess molar entropy [J/K/mol]. Property specific to solutions.

property excess_gibbs

Returns molar excess gibbs free energy [J/mol]. Property specific to solutions.

property excess_partial_entropies

Returns excess partial entropies [J/K]. Property specific to solutions.

property excess_partial_gibbs

Returns excess partial molar gibbs free energy [J/mol]. Property specific to solutions.

property excess_partial_volumes

Returns excess partial volumes [m^3]. Property specific to solutions.

property excess_volume

Returns excess molar volume of the solution [m^3/mol]. Specific property for solutions.

property formula

Returns molar chemical formula of the solution. :rtype: Counter

property gibbs

Alias for molar_gibbs()

property gibbs_hessian

Returns an array containing the second compositional derivative of the Gibbs free energy [J]. Property specific to solutions.

property gr

Alias for grueneisen_parameter()

property grueneisen_parameter

Returns grueneisen parameter of the solution [unitless]. Aliased with self.gr.

property helmholtz

Alias for molar_helmholtz()

property independent_element_indices

A list of an independent set of element indices. If the amounts of these elements are known (element_amounts), the amounts of the other elements can be inferred by -compositional_null_basis[independent_element_indices].dot(element_amounts).

property isentropic_bulk_modulus_reuss

Returns adiabatic bulk modulus of the solution [Pa]. Aliased with self.K_S.

property isentropic_compressibility_reuss

Returns adiabatic compressibility of the solution. (or inverse adiabatic bulk modulus) [1/Pa]. Aliased with self.K_S.

property isentropic_thermal_gradient

Returns:

dTdP, the change in temperature with pressure at constant entropy [Pa/K]

Return type:

float

property isothermal_bulk_modulus_reuss

Returns isothermal bulk modulus of the solution [Pa]. Aliased with self.K_T.

property isothermal_compressibility_reuss

Returns isothermal compressibility of the solution. (or inverse isothermal bulk modulus) [1/Pa]. Aliased with self.K_T.

property molar_enthalpy

Returns molar enthalpy of the solution [J/mol]. Aliased with self.H.

property molar_entropy

Returns molar entropy of the solution [J/K/mol]. Aliased with self.S.

property molar_gibbs

Returns molar Gibbs free energy of the solution [J/mol]. Aliased with self.gibbs.

property molar_heat_capacity_p

Returns molar heat capacity at constant pressure of the solution [J/K/mol]. Aliased with self.C_p.

property molar_heat_capacity_v

Returns molar heat capacity at constant volume of the solution [J/K/mol]. Aliased with self.C_v.

property molar_helmholtz

Returns molar Helmholtz free energy of the solution [J/mol]. Aliased with self.helmholtz.

property molar_internal_energy

Returns molar internal energy of the mineral [J/mol]. Aliased with self.energy

property molar_mass

Returns molar mass of the solution [kg/mol].

property molar_volume

Returns molar volume of the solution [m^3/mol]. Aliased with self.V.

property n_endmembers

The number of endmembers in the solution.

property n_reactions

The number of reactions in reaction_basis.

property name

Human-readable name of this material.

By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.

property p_wave_velocity

Returns P wave speed of the solution [m/s]. Aliased with self.v_p.

property partial_entropies

Returns endmember partial entropies [J/K]. Property specific to solutions.

property partial_gibbs

Returns endmember partial molar gibbs free energy [J/mol]. Property specific to solutions.

property partial_volumes

Returns endmember partial volumes [m^3]. Property specific to solutions.

property pressure

Returns current pressure that was set with set_state().

Note

Aliased with P().

Returns:

Pressure in [Pa].

Return type:

float

print_minerals_of_current_state()

Print a human-readable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.

property reaction_basis

An array where each element arr[i,j] corresponds to the number of moles of endmember[j] involved in reaction[i].

reset()

Resets all cached material properties.

It is typically not required for the user to call this function.

property rho

Alias for density()

set_composition(molar_fractions)

Set the composition for this solution. Resets cached properties.

Parameters:

molar_fractions (list of float) – Molar abundance for each endmember, needs to sum to one.

set_method(method)

Set the equation of state to be used for this mineral. Takes a string corresponding to any of the predefined equations of state: ‘bm2’, ‘bm3’, ‘mgd2’, ‘mgd3’, ‘slb2’, ‘slb3’, ‘mt’, ‘hp_tmt’, or ‘cork’. Alternatively, you can pass a user defined class which derives from the equation_of_state base class. After calling set_method(), any existing derived properties (e.g., elastic parameters or thermodynamic potentials) will be out of date, so set_state() will need to be called again.

set_state(pressure, temperature)

(copied from set_state):

Set the material to the given pressure and temperature.

Parameters:

• pressure (float) – The desired pressure in [Pa].

• temperature (float) – The desired temperature in [K].

set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])

This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure self-consistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq root-finding method.

Parameters:

• volume (float) – The desired molar volume of the mineral [m^3].

• temperature (float) – The desired temperature of the mineral [K].

• pressure_guesses (list) – A list of floats denoting the initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state. Defaults to [0.e9, 10.e9].

property shear_modulus

Returns shear modulus of the solution [Pa]. Aliased with self.G.

property shear_wave_velocity

Returns shear wave speed of the solution [m/s]. Aliased with self.v_s.

site_formula(precision=2)

Returns the molar chemical formula of the solution with

site occupancies. For example, [Mg0.4Fe0.6]2SiO4.

Parameters:

precision (int) – Precision with which to print the site occupancies

Returns:

Molar chemical formula of the solution with site occupancies

Return type:

str

property site_occupancies

Returns:

The fractional occupancies of species on each site.

Return type:

list of OrderedDicts

property stoichiometric_array

An array where each element arr[i,j] corresponds to the number of atoms of element[j] in endmember[i].

property stoichiometric_matrix

A sympy Matrix where each element M[i,j] corresponds to the number of atoms of element[j] in endmember[i].

property temperature

Returns current temperature that was set with set_state().

Note

Aliased with T().

Returns:

Temperature in [K].

Return type:

float

property thermal_expansivity

Returns thermal expansion coefficient (alpha) of the solution [1/K]. Aliased with self.alpha.

to_string()

Returns the name of the mineral class

unroll()

Unroll this material into a list of burnman.Mineral and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure).

Note

Needs to be implemented in derived classes.

Returns:

A list of molar fractions which should sum to 1.0, and a list of burnman.Mineral objects containing the minerals in the material.

Return type:

tuple

property v_p

Alias for p_wave_velocity()

property v_phi

Alias for bulk_sound_velocity()

property v_s

Alias for shear_wave_velocity()

property volume_hessian

Returns an array containing the second compositional derivative of the volume [m^3]. Property specific to solutions.

burnman.optimize.composition_fitting.fit_composition_to_solution(solution, fitted_variables, variable_values, variable_covariances, variable_conversions=None, normalize=True)

Takes a Solution object and a set of variable names and associates values and covariances and finds the molar fractions of the solution which provide the best fit (in a least-squares sense) to the variable values.

The fitting applies appropriate non-negativity constraints (i.e. no species can have a negative occupancy on a site).

Parameters:

• solution (burnman.Solution) – The solution to use in the fitting procedure.

• fitted_variables (list of str) – A list of the variables used to find the best-fit molar fractions of the solution. These should either be elements such as “Fe”, site_species such as “Fef_B” which would correspond to a species labelled Fef on the second site, or user-defined variables which are arithmetic sums of elements and/or site_species defined in “variable_conversions”.

• variable_values (numpy.array) – Numerical values of the fitted variables. These should be given as amounts; they do not need to be normalized.

• variable_covariances (2D numpy.array) – Covariance matrix of the variables.

• variable_conversions (dict of dict, or None) – A dictionary converting any user-defined variables into an arithmetic sum of element and site-species amounts. For example, {‘Mg_equal’: {‘Mg_A’: 1., ‘Mg_B’: -1.}}, coupled with Mg_equal = 0 would impose a constraint that the amount of Mg would be equal on the first and second site in the solution.

• normalize (bool) – If True, normalizes the optimized molar fractions to sum to unity.

Returns:

Optimized molar fractions, corresponding covariance matrix and the weighted residual.

Return type:

tuple of 1D numpy.array, 2D numpy.array and float

burnman.optimize.composition_fitting.fit_phase_proportions_to_bulk_composition(phase_compositions, bulk_composition)

Performs weighted constrained least squares on a set of phase compositions to find the amount of those phases that best-fits a given bulk composition.

The fitting applies appropriate non-negativity constraints (i.e. no phase can have a negative abundance in the bulk).

Parameters:

• phase_compositions (2D numpy.array) – The composition of each phase. Can be in weight or mole amounts.

• bulk_composition (numpy.array) – The bulk composition of the composite. Must be in the same units as the phase compositions.

Returns:

Optimized molar fractions, corresponding covariance matrix and the weighted residual.

Return type:

tuple of 1D numpy.array, 2D numpy.array and float

class burnman.optimize.eos_fitting.MineralFit(mineral, data, data_covariances, flags, fit_params, mle_tolerances, delta_params=None, bounds=None)

Bases: object

Class for fitting mineral parameters to experimental data. Instances of this class are passed to burnman.nonlinear_least_squares_fit().

For attributes added to this model when fitting is done, please see the documentation for that function.

set_params(param_values)

get_params()

function(x, flag)

normal(x, flag)

burnman.optimize.eos_fitting.fit_PTp_data(mineral, fit_params, flags, data, data_covariances=[], mle_tolerances=[], param_tolerance=1e-05, delta_params=None, bounds=None, max_lm_iterations=50, verbose=True)

Given a mineral of any type, a list of fit parameters and a set of P-T-property points and (optional) uncertainties, this function returns a list of optimized parameters and their associated covariances, fitted using the scipy.optimize.curve_fit routine.

Parameters:

• mineral (burnman.Mineral) – Mineral for which the parameters should be optimized.

• fit_params (list of str) – List of dictionary keys contained in mineral.params corresponding to the variables to be optimized during fitting. Initial guesses are taken from the existing values for the parameters

• flags (string or list of strings) – Attribute names for the property to be fit for the whole dataset or each datum individually (e.g. ‘V’)

• data (2D numpy.array) – Observed X-P-T-property values

• data_covariances (3D numpy.array) – X-P-T-property covariances (optional) If not given, all covariance matrices are chosen such that all data points have equal weight, with all error in the pressure.

• mle_tolerances (numpy.array) – Tolerances for termination of the maximum likelihood iterations (optional).

• param_tolerance (float) – Fractional tolerance for termination of the nonlinear optimization (optional).

• delta_params (numpy.array) – Initial values for the change in parameters (optional).

• bounds (2D numpy.array) – Minimum and maximum bounds for the parameters (optional). The shape must be (n_parameters, 2).

• max_lm_iterations (int) – Maximum number of Levenberg-Marquardt iterations.

• verbose (bool) – Whether to print detailed information about the optimization to screen.

Returns:

Model with optimized parameters.

Return type:

burnman.optimize.eos_fitting.MineralFit

burnman.optimize.eos_fitting.fit_PTV_data(mineral, fit_params, data, data_covariances=[], delta_params=None, bounds=None, param_tolerance=1e-05, max_lm_iterations=50, verbose=True)

A simple alias for the fit_PTp_data for when all the data is volume data

class burnman.optimize.eos_fitting.SolutionFit(solution, data, data_covariances, flags, fit_params, mle_tolerances, delta_params=None, bounds=None)

Bases: object

Class for fitting mineral parameters to experimental data. Instances of this class are passed to burnman.nonlinear_least_squares_fit().

For attributes added to this model when fitting is done, please see the documentation for that function.

set_params(param_values)

get_params()

function(x, flag)

normal(x, flag)

burnman.optimize.eos_fitting.fit_XPTp_data(solution, fit_params, flags, data, data_covariances=[], mle_tolerances=[], param_tolerance=1e-05, delta_params=None, bounds=None, max_lm_iterations=50, verbose=True)

Given a symmetric solution, a list of fit parameters and a set of P-T-property points and (optional) uncertainties, this function returns a list of optimized parameters and their associated covariances, fitted using the scipy.optimize.curve_fit routine.

Parameters:

• solution (burnman.Solution) – Solution for which the parameters should be optimized.

• fit_params (list of lists) – Variables to be optimized during fitting. Each list is either of length two or three. The first item of length-2 lists should be a dictionary key contained in one of the endmember mineral.params, and the second item should be the index of the endmember in the solution (indexing starts from 0). The first item of length-3 lists should be one of ‘E’, ‘S’ or ‘V’ (the excess energies, entropies or volumes in each binary). The second two items should be the indices of the pair of endmembers bounding the binary, in ascending order (indexing starts from 0). Initial guesses are taken from the existing values for the parameters.

• flags (string or list of strings) – Attribute names for the property to be fit for the whole dataset or each datum individually (e.g. ‘V’)

• data (2D numpy.array) – Observed X-P-T-property values

• data_covariances (3D numpy.array) – X-P-T-property covariances (optional). If not given, all covariance matrices are chosen such that all data points have equal weight, with all error in the pressure.

• mle_tolerances (numpy.array) – Tolerances for termination of the maximum likelihood iterations (optional).

• param_tolerance (float) – Fractional tolerance for termination of the nonlinear optimization (optional).

• delta_params (numpy.array) – Initial values for the change in parameters (optional).

• bounds (2D numpy.array) – Minimum and maximum bounds for the parameters (optional). The shape must be (n_parameters, 2).

• max_lm_iterations (int) – Maximum number of Levenberg-Marquardt iterations.

• verbose (bool) – Whether to print detailed information about the optimization to screen.

Returns:

Model with optimized parameters.

Return type:

burnman.optimize.eos_fitting.SolutionFit

burnman.optimize.linear_fitting.weighted_constrained_least_squares(A, b, Cov_b=None, equality_constraints=None, inequality_constraints=None)

Solves a weighted, constrained least squares problem using cvxpy. The objective function is to minimize the following: sum_squares(Cov_b^(-1/2).A.x - Cov_b^(-1/2).b)) subject to C.x == c D.x <= d

Parameters:

• A (2D numpy array) – An array defining matrix A in the objective function above.

• b (numpy array) – An array defining vector b in the objective function above.

• Cov_b (2D numpy array) – A covariance matrix associated with b

• equality_constraints (list containing a 2D array and 1D array) – A list containing the matrices C and c in the objective function above.

• inequality_constraints (list containing a 2D array and 1D array) – A list containing the matrices D and d in the objective function above.

Returns:

Tuple containing the optimized phase amounts (1D numpy.array), a covariance matrix corresponding to the optimized phase amounts (2D numpy.array), and the weighted residual of the fitting procedure (a float).

Return type:

tuple

burnman.optimize.nonlinear_fitting.nonlinear_least_squares_fit(model, lm_damping=0.0, param_tolerance=1e-07, max_lm_iterations=100, verbose=False)

Function to compute the “best-fit” parameters for a model by nonlinear least squares fitting.

The nonlinear least squares algorithm closely follows the logic in Section 23.1 of Bayesian Probability Theory (von der Linden et al., 2014; Cambridge University Press).

Parameters:

:param model: Model containing data to be fit, and functions to

aid in fitting.

:type model: object

:param lm_damping: Levenberg-Marquardt parameter for least squares minimization.

:type lm_damping: float

:param param_tolerance: Levenberg-Marquardt iterations are terminated when

the maximum fractional change in any of the parameters during an iteration drops below this value

:type param_tolerance: float

:param max_lm_iterations: Maximum number of Levenberg-Marquardt iterations

:type max_lm_iterations: int

:param verbose: Print some information to standard output

:type verbose: bool

.. note:: The object passed as model must have the following attributes:

• data [2D numpy.array] - Elements of x[i][j] contain the

  observed position of data point i.

• data_covariances [3D numpy.array] Elements of cov[i][j][k] contain

  the covariance matrix of data point i.

• mle_tolerances [numpy.array] - The iterations to find the maximum likelihood

  estimator for each observed data point will stop when mle_tolerances[i] < np.linalg.norm(data_mle[i] - model.function(data_mle[i], flag))

• delta_params [numpy.array] - parameter perturbations used to compute the jacobian

Must also have the following methods: * set_params(self, param_values) - Function to set parameters. * get_params(self) - Function to get current model parameters. * function(self, x) - Returns value of model function evaluated at x. * normal(self, x) - Returns value of normal to the model function evaluated at x.

After this function has been performed, the following attributes are added to model:

• n_dof [int] - Degrees of freedom of the system.

• data_mle [2D numpy array] - Maximum likelihood estimates of the observed data points

  on the best-fit curve.

• jacobian [2D numpy array] - d(weighted_residuals)/d(parameter).

• weighted_residuals [numpy array] - Weighted residuals.

• weights [numpy array] - 1/(data variances normal to the best fit curve).

• WSS [float] - Weighted sum of squares residuals.

• popt [numpy array] - Optimized parameters.

• pcov [2D numpy array] - Covariance matrix of optimized parameters.

• noise_variance [float] - Estimate of the variance of the data normal to the curve.

This function is available as ``burnman.nonlinear_least_squares_fit``.

burnman.optimize.nonlinear_fitting.confidence_prediction_bands(model, x_array, confidence_interval, f, flag=None)

This function calculates the confidence and prediction bands of the function f(x) from a best-fit model with uncertainties in its parameters as calculated (for example) by the function nonlinear_least_squares_fit().

The values are calculated via the delta method, which estimates the variance of f evaluated at x as var(f(x)) = df(x)/dB var(B) df(x)/dB where df(x)/dB is the vector of partial derivatives of f(x) with respect to B.

Parameters:

• model (object) – As modified (for example) by the function burnman.nonlinear_least_squares_fit(). Should contain the following functions: get_params, set_params, function, normal And attributes: delta_params, pcov, dof, noise_variance

• x_array (2D numpy.array) – Coordinates at which to evaluate the bounds.

• confidence_interval (float) – Probability level of finding the true model (confidence bound) or any new data point (probability bound). For example, the 95% confidence bounds should be calculated using a confidence interval of 0.95.

• f (function) – The function defining the variable y=f(x) for which the confidence and prediction bounds are desired.

• flag (type informed by model object) – This (optional) flag is passed to model.function to control how the modified position of x is calculated. This value is then used by f(x)

Returns:

An element of bounds[i][j] gives the lower and upper confidence (i=0, i=1) and prediction (i=2, i=3) bounds for the jth data point.

Return type:

2D numpy.array

burnman.optimize.nonlinear_fitting.abs_line_project(M, n)

burnman.optimize.nonlinear_fitting.plot_cov_ellipse(cov, pos, nstd=2, ax=None, **kwargs)

Plots an nstd sigma error ellipse based on the specified covariance matrix (cov). Additional keyword arguments are passed on to the ellipse patch artist.

Parameters:

• cov (numpy.array) – The 2x2 covariance matrix to base the ellipse on.

• pos (list or numpy.array) – The location of the center of the ellipse. Expects a 2-element sequence of [x0, y0].

• nstd (float) – The radius of the ellipse in numbers of standard deviations. Defaults to 2 standard deviations.

• ax (matplotlib.pyplot.axes) – The axis that the ellipse will be plotted on. Defaults to the current axis.

• kwargs – Additional keyword arguments are passed on to the ellipse patch.

Returns:

The covariance ellipse (already applied to the desired axes object).

Return type:

matplotlib.patches.Ellipse

burnman.optimize.nonlinear_fitting.corner_plot(popt, pcov, param_names=[], n_std=1.0)

Creates a corner plot of covariances

Parameters:

• popt (numpy.array) – Optimized parameters.

• pcov (2D numpy.array) – Covariance matrix of the parameters.

• param_names (list) – Parameter names.

• n_std (float) – Number of standard deviations for ellipse.

Returns:

matplotlib.pyplot.figure and list of matplotlib.pyplot.Axes objects.

Return type:

tuple

burnman.optimize.nonlinear_fitting.weighted_residual_plot(ax, model, flag=None, sd_limit=3, cmap=<matplotlib.colors.LinearSegmentedColormap object>, plot_axes=[0, 1], scale_axes=[1.0, 1.0])

Creates a plot of the weighted residuals The user can choose the projection axes, and scaling to apply to those axes The chosen color palette (cmap) is discretised by standard deviation up to a cut off value of sd_limit.

Parameters:

• ax – Plot.

• type – matplotlib.pyplot.Axes

• model (object) – A model as used by burnman.nonlinear_least_squares_fit(). Must contain the attributes model.data, model.weighted_residuals and model.flags (if flag is not None).

• flag (str) – String to determine which data to plot. Finds matches with model.flags.

• sd_limit (float) – Data with weighted residuals exceeding this limit are plotted in black.

• cmap (matplotlib color palette) – Color palette.

• plot_axes (list of int) – Data axes to use as plot axes.

• scale_axes (list of float) – Plot axes are scaled by multiplication of the data by these values.

Returns:

Coloured scatter plot of the weighted residuals in data space.

Return type:

matplotlib Axes object

burnman.optimize.nonlinear_fitting.extreme_values(weighted_residuals, confidence_interval)

This function uses extreme value theory to calculate the number of standard deviations away from the mean at which we should expect to bracket all of our n data points at a certain confidence level.

It then uses that value to identify which (if any) of the data points lie outside that region, and calculates the corresponding probabilities of finding a data point at least that many standard deviations away.

Parameters:

• weighted_residuals (array of float) – Array of residuals weighted by the square root of their variances wr_i = r_i/sqrt(var_i).

• confidence_interval (float) – Probability at which all the weighted residuals lie within the confidence bounds.

Returns:

Number of standard deviations at which we should expect to encompass all data at the user-defined confidence interval, the indices of weighted residuals exceeding the confidence_interval defined by the user, and the probabilities that the extreme data point of the distribution lies further from the mean than the observed position wr_i for each i in the “indices” output array.

Return type:

tuple of (float, numpy.array, numpy.array)

burnman.optimize.nonlinear_fitting.plot_residuals(ax, weighted_residuals, n_bins=None, flags=[])

burnman.optimize.nonlinear_solvers.solve_constraint_lagrangian(x, jac_x, c_x, c_prime)

Function which solves the problem minimize || J.dot(x_mod - x) || subject to C(x_mod) = 0 via the method of Lagrange multipliers.

Parameters:

• x (1D numpy array) – Parameter values at x.

• jac_x (2D numpy array.) – The (estimated, approximate or exact) value of the Jacobian J(x).

• c_x (1D numpy array) – Values of the constraints at x.

• c_prime (2D array of floats) – The Jacobian of the constraints (A, where A.x + b = 0).

Returns:

An array containing the parameter values which minimize the L2-norm of any function which has the Jacobian jac_x, and another array containing the multipliers for each of the equality constraints.

Return type:

tuple

burnman.optimize.nonlinear_solvers.damped_newton_solve(F, J, guess, tol=1e-06, max_iterations=100, lambda_bounds=<function <lambda>>, linear_constraints=(0.0, array([-1.])), store_iterates=False)

Solver for the multivariate nonlinear system F(x)=0 with Jacobian J(x), using the damped affine invariant modification to Newton’s method (Deuflhard, 1974;1975;2004). Here we follow the algorithm as described in Nowak and Weimann (1991): [Technical Report TR-91-10, Algorithm B], modified to accept linear inequality constraints.

Linear inequality constraints are provided by the arrays constraints_A and constraints_b. The constraints are satisfied if A*x + b <= 0. If any constraints are not satisfied by the current value of lambda, lambda is reduced to satisfy all the constraints.

If a current iterate starting point (x_i) lies on one or more constraints and the Newton step violates one or more of those constraints, then the next step is calculated via the method of Lagrangian multipliers, minimizing the L2-norm of F(x_i+1) subject to the violated constraints.

Successful termination of the solver is based on three criteria:

• all(np.abs(dx (simplified newton step) < tol))

• all(np.abs(dx (full Newton step) < sqrt(10*tol))) (avoids pathology) and

• lambda = lambda_bounds(dx, x)[1] (lambda = 1 for a full Newton step).

If these criteria are not satisfied, iterations continue until one of the following occurs:

• the value of lmda is reduced to its minimum value (for v. nonlinear problems)

• successive iterations have descent vectors which violate the constraints

• the maximum number of iterations (given by max_iterations) is reached.

Information on the root (or lack of root) obtained by the solver is provided in the returned namedtuple.

Parameters:

• F (function of x) – Function returning the system function F(x) as a 1D numpy array.

• J (function of x) – Function returning the Jacobian function J(x) as a 2D numpy array.

• guess (1D numpy.array) – Starting guess for the solver.

• tol (float or array of floats) – Tolerance(s) for termination.

• max_iterations (int) – Maximum number of iterations for the solver.

• lambda_bounds – A function of dx and x that returns a tuple of floats corresponding to the minimum and maximum allowed fractions of the full newton step (dx).

• linear_constraints (tuple of a 2D numpy.array (A) and 1D numpy.array (b)) – tuple of a 2D numpy array (A) and 1D numpy array (b) Constraints are satisfied if A.x + b < eps

Returns:

A namedtuple with the following attributes:

• x: The solution vector [1D numpy array].

• F: The evaluated function F(x) [1D numpy array].

• F_norm: Euclidean norm of F(x) [float].

• J: The evaluated Jacobian J(x) [2D numpy array].

• n_it: Number of iterations [int].

• code: Numerical description of the solver termination [int].

  • 0: Successful convergence

  • 1: Failure due to solver hitting lower lambda bound

  • 2: Failure due to descent vector crossing constraints

  • 3: Failure due to solver reaching maximum number of iterations

• text: Description of the solver termination [str].

• success: Solver convergence state [bool].

• iterates: [namedtuple]

  Only present if store_iterates=True Includes the following attributes:

  • x: list of 1D numpy arrays of floats

    The parameters for each iteration

  • F: list of 2D numpy arrays of floats

    The function for each iteration

  • lmda: list of floats

    The value of the damping parameter for each iteration

Return type:

namedtuple