Optimization¶

class burnman.optimize.composition_fitting.DummyCompositionSolution(endmember_element_formulae, endmember_site_formulae)[source]¶

Bases: Solution

A mock Solution class to enable compositional fitting when no thermodynamic model is available. Useful for pure composition analysis.

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 endmembers in the solution.

property C_p¶: Alias for heat_capacity_p()

property C_v¶: Alias for heat_capacity_v()

property G¶: Alias for shear_modulus()

property G_eff¶: Alias for effective_shear_modulus()

property H¶: Alias for enthalpy()

property K_S¶: Alias for isentropic_bulk_modulus_reuss()

property K_T¶: Alias for isothermal_bulk_modulus_reuss()

property K_eff¶: Alias for effective_isentropic_bulk_modulus()

property P¶: Alias for pressure()

property S¶: Alias for entropy()

property T¶: Alias for temperature()

property V¶: Alias for 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¶

Returns:

A basis for the compositional degrees of freedom, orthogonal to the reaction space.

Each row of the returned array defines an independent compositional constraint in terms of moles of endmembers. These constraints span the row space of the stoichiometric matrix and are complementary to the reaction basis (which spans the left null space).

For example, endmembers with site-species occupancies [Fe][Mg] and [Mg][Fe] have reaction basis vector [1, -1]. A compositional basis vector orthogonal to this is: [1, 1], which satisfies the requirement that the Fe-Mg ratio of the system remains fixed and the total amount of the two sites remains equal.

Return type:

numpy.ndarray

property compositional_null_basis¶

Returns a basis for the compositional nullspace (kernel) of the stoichiometric matrix.

This basis consists of vectors that represent linear constraints any valid bulk composition must satisfy to be formed from the endmembers. In other words, for any bulk composition vector b that can be expressed as a combination of the endmembers, multiplying it by these vectors results in zero (N @ b = 0).

These constraints typically reflect conserved properties such as elemental mass balance or charge neutrality.

Example: Suppose the stoichiometric matrix corresponds to elements C and O in species, and the nullspace basis vector is [1, -2]. Then for any bulk composition vector b = [b_C, b_O], the constraint 1 * b_C - 2 * b_O = 0 must hold. This means the ratio of carbon to oxygen must satisfy that relationship to be achievable from the given endmembers.

Returns:: A 2D array where each row is a basis vector of the nullspace.
Return type:: numpy.ndarray

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 effective_isentropic_bulk_modulus¶

Returns the effective isentropic bulk modulus of the mineral.

Note

Needs to be implemented in derived classes. Aliased with K_eff().

Returns:: Effective isentropic bulk modulus in [Pa].
Return type:: float

property effective_shear_modulus¶

Returns the effective shear modulus of the mineral.

Note

Needs to be implemented in derived classes. Aliased with G_eff().

Returns:: Effective shear modulus in [Pa].
Return type:: float

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 endmembers in the solution.

property endmembers¶

property enthalpy¶: Extensive enthalpy of the material (J).

property entropy¶: Extensive entropy of the material (J/K).

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:

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¶: Extensive Gibbs free energy of the material (J).

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 heat_capacity_p¶: Extensive heat capacity at constant pressure of the material (J/K).

property heat_capacity_v¶: Extensive heat capacity at constant volume of the material (J/K).

property helmholtz¶: Extensive Helmholtz free energy of the material (J).

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 internal_energy¶: Extensive internal energy of the material (J).

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 mass¶

Returns the mass of this material.

Returns:: Mass in [kg].
Return type:: float

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 number_of_moles¶

Returns the number of moles in this material.

Returns:: Number of moles.
Return type:: float

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¶

Returns a basis for the reaction space, represented as the left nullspace of the stoichiometric matrix.

Each row of the returned array corresponds to a linearly independent chemical reaction that conserves mass and charge. The value at arr[i, j] gives the number of moles of endmember[j] involved in reaction[i].

Returns:: An array where each row defines a balanced chemical reaction in terms of molar amounts of endmembers.
Return type:: numpy.ndarray

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: ‘bm3shear2’, ‘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.

If the mineral is an endmember that has a pressure(T, V) function and does not have property modifiers, this function evaluates the pressure directly. Otherwise, it finds the pressure using the brentq root-finding method on the molar_volume attribute of the mineral. This ensures self-consistency even if the mineral has a pressure-dependent volume modifier.

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, print_absent_species=True, site_occupancies=None)¶

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
print_absent_species (bool) – If True, prints site occupancies that are absent (i.e., zero).
site_occupancies (numpy.ndarray) – Optional site occupancies to use instead of the current ones.

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 as a flat vector.
Return type:: numpy.ndarray

property stoichiometric_array¶: A numpy 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¶: Extensive volume of the material (m^3).

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)[source]¶

Fit a set of molar fractions to a solution model that best matches observed compositional variables (in a least-squares sense).

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, default=True) – Whether to normalize the molar fractions to sum to one.

Returns:

(molar_fractions, covariance_matrix, weighted_residual):

molar_fractions: Optimized molar fractions of each endmember.
covariance_matrix: Covariance matrix of the optimized molar
fractions, computed from the pseudo-inverse of the weighted design matrix.
weighted_residual: The weighted residual of the fit, representing
the square root of the weighted sum of squared differences between observed and modeled variables.

Return type:

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

See also

burnman.optimize.linear_fitting.weighted_constrained_least_squares() - Performs the weighted constrained least squares optimization used internally.

burnman.optimize.eos_fitting.default_mle_tolerances(material, flags)[source]¶

Computes default tolerances for maximum likelihood estimation (MLE).

The tolerances are computed per property flag, using either: - A fixed absolute tolerance (e.g., 1.0 J) for energy-related properties - A fixed absolute tolerance (e.g., 100 Pa) for pressure - A relative tolerance (1e-5 x standard state value) for all other properties

All computed tolerances must be > 0.

Parameters:

material (burnman.Material) – The material instance to evaluate properties on.
flags (list of str) – List of property names (strings) for which MLE tolerances should be generated.

Returns:

A NumPy array of tolerances, one per property flag.

Return type:

np.ndarray

Raises:

Exception – If any computed tolerance is not strictly positive.

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

Bases: NonLinearModel

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)[source]¶

Parameters:: param_values (numpy.ndarray) – Parameter values to set.

get_params()[source]¶

Returns:: Current model parameters as a NumPy array.
Return type:: numpy.ndarray

function(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Model output.

Return type:

numpy.ndarray

normal(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Normal vector.

Return type:

numpy.ndarray

validate()¶

Ensure that all required model attributes are present before fitting. These include observed data, data covariance matrices, convergence tolerances for MLE projections, and finite difference step sizes.

Raises:: AssertionError – If required attributes are missing.

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, param_priors=None, param_prior_inv_cov_matrix=None, verbose=True)[source]¶

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, param_priors=None, param_prior_inv_cov_matrix=None, verbose=True)[source]¶: A simple alias for the fit_PTp_data for when all the data is volume data

class burnman.optimize.eos_fitting.MineralFitV(mineral, data, data_covariances, flags, fit_params, mle_tolerances, delta_params=None, bounds=None)[source]¶

Bases: NonLinearModel

Class for fitting mineral parameters to experimental data using volume as an independent variable. 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)[source]¶

Parameters:: param_values (numpy.ndarray) – Parameter values to set.

get_params()[source]¶

Returns:: Current model parameters as a NumPy array.
Return type:: numpy.ndarray

function(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Model output.

Return type:

numpy.ndarray

normal(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Normal vector.

Return type:

numpy.ndarray

validate()¶

Raises:: AssertionError – If required attributes are missing.

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

Given a mineral of any type, a list of fit parameters and a set of V-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. ‘P’)
data (2D numpy.array) – Observed V-T-property values
data_covariances (3D numpy.array) – V-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 properties.
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.MineralFitV

burnman.optimize.eos_fitting.fit_VTP_data(mineral, fit_params, data, data_covariances=[], delta_params=None, bounds=None, param_tolerance=1e-05, max_lm_iterations=50, param_priors=None, param_prior_inv_cov_matrix=None, verbose=True)[source]¶: A simple alias for the fit_VTp_data for when all the data is pressure data

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

Bases: NonLinearModel

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)[source]¶

Parameters:: param_values (numpy.ndarray) – Parameter values to set.

get_params()[source]¶

Returns:: Current model parameters as a NumPy array.
Return type:: numpy.ndarray

function(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Model output.

Return type:

numpy.ndarray

validate()¶

Raises:: AssertionError – If required attributes are missing.

normal(x, flag)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Normal vector.

Return type:

numpy.ndarray

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, param_priors=None, param_prior_inv_cov_matrix=None, verbose=True)[source]¶

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, allow_rank_deficient=False)[source]¶

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.
allow_rank_deficient (bool) – If True, allows the problem to be solved even if the design matrix is rank-deficient.

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

class burnman.optimize.nonlinear_fitting.NonLinearModel[source]¶

Bases: ABC

Abstract base class that defines the required interface for models used in nonlinear least squares fitting. Models must support evaluation, normal vector computation, and parameter handling.

abstractmethod get_params()[source]¶

Returns:: Current model parameters as a NumPy array.
Return type:: numpy.ndarray

abstractmethod set_params(param_values)[source]¶

Parameters:: param_values (numpy.ndarray) – Parameter values to set.

abstractmethod function(x, flag=None)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Model output.

Return type:

numpy.ndarray

abstractmethod normal(x, flag=None)[source]¶

Parameters:

x (numpy.ndarray) – Input coordinates.
flag (any) – Optional argument to control evaluation mode.

Returns:

Normal vector.

Return type:

numpy.ndarray

validate()[source]¶

Raises:: AssertionError – If required attributes are missing.

burnman.optimize.nonlinear_fitting.abs_line_project(M, n)[source]¶

Project a covariance matrix M onto a unit direction vector n. This gives the variance of the distribution in the direction n.

Parameters:

M (numpy.ndarray) – Covariance matrix.
n (numpy.ndarray) – Direction vector.

Returns:

Projected variance.

Return type:

float

burnman.optimize.nonlinear_fitting.mle_estimate(model, x, x_m, cov, flag)[source]¶

Find an approximation to the maximum likelihood estimate of a point on the model surface corresponding to noisy observation x, under Gaussian noise with covariance cov, given a starting estimate on the model surface x_m.

This approximation is the orthogonal projection of x onto the local tangent plane of the model surface, as determined using the Mahalanobis metric defined by cov.

Parameters:

model (NonLinearModel) – Model object.
x (numpy.ndarray) – Observed data point.
x_m (numpy.ndarray) – Initial guess on the model surface.
cov (numpy.ndarray) – Covariance matrix.
flag (any) – Optional flag for evaluation.

Returns:

Tuple containing MLE position, residual, and projected variance.

Return type:

tuple(numpy.ndarray, float, float)

burnman.optimize.nonlinear_fitting.find_mle(model)[source]¶

Find the maximum likelihood point for each datum given the current model parameters.

This is done iteratively using the function mle_estimate.

Parameters:: model (FittableModel) – Model object with required attributes.
Returns:: MLE data positions, weighted residuals, and weights.
Return type:: tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray)

burnman.optimize.nonlinear_fitting.calculate_jacobian(model)[source]¶

Compute the Jacobian matrix of weighted residuals with respect to model parameters using central finite differences.

Parameters:: model (FittableModel) – Model object.
Modifies:: model.jacobian — Populated with partial derivatives.

burnman.optimize.nonlinear_fitting.nonlinear_least_squares_fit(model, lm_damping=0.0, param_tolerance=1e-07, max_lm_iterations=100, param_priors=None, param_prior_inv_cov_matrix=None, verbose=False)[source]¶

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:

model (FittableModel) – Model with fitting interface.
lm_damping (float) – Damping factor for Levenberg-Marquardt updates.
param_tolerance (float) – Convergence tolerance based on fractional parameter change.
max_lm_iterations (int) – Maximum number of LM iterations.
param_priors (1D numpy array) – Prior values for the parameters.
param_prior_inv_cov_matrix (2D numpy array (square)) – Inverse of the 1 sigma uncertainties for the prior values of the parameters.
verbose (bool) – If True, print iteration status.

Modifies:

model — Sets optimized parameters, covariance matrix, weighted residuals, Jacobian, and noise estimates.

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)[source]¶

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.plot_cov_ellipse(cov, pos, nstd=2, ax=None, **kwargs)[source]¶

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)[source]¶

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])[source]¶

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)[source]¶

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=[])[source]¶

class burnman.optimize.nonlinear_solvers.Iterates(x: list = <factory>, F: list = <factory>, lmda: list = <factory>)[source]¶

Bases: object

Container for storing iteration history of the solver. Attributes are lists of the iterates x, function evaluations F, and step lengths relative to the full Newton step lmda.

x: list¶

F: list¶

lmda: list¶

append(x, F, lmda)[source]¶: Append an iteration’s data to the history. :param x: Current iterate. :param F: Current function evaluation. :param lmda: Current step length.

close()[source]¶: Finalize the iteration history by converting lists to arrays.

class burnman.optimize.nonlinear_solvers.TerminationCode(value)[source]¶

Bases: Enum

Enumeration of termination codes for the nonlinear solver. These are used to indicate the reason for solver termination. The codes are:

SUCCESS (0): The solver successfully found a root.

TOO_NONLINEAR (1): The function is too non-linear for lower lambda bound.

CONSTRAINT_VIOLATION (2): The descent vector crosses the constraints.

MAX_ITERATIONS (3): The solver reached the maximum number of iterations.

SINGULAR_FAIL (4): The system was effectively singular and failed to converge.

SINGULAR_SUCCESS (5): The solver found a root but the system is singular.

SUCCESS = 0¶

TOO_NONLINEAR = 1¶

CONSTRAINT_VIOLATION = 2¶

MAX_ITERATIONS = 3¶

SINGULAR_FAIL = 4¶

SINGULAR_SUCCESS = 5¶

class burnman.optimize.nonlinear_solvers.Solution(x: ndarray | None = None, n_it: int = 0, F: ndarray | None = None, F_norm: float | None = None, J: ndarray | None = None, code: TerminationCode | None = None, text: str | None = None, success: bool = False, store_iterates: bool = False)[source]¶

Bases: object

Container for the solution of a nonlinear system of equations.

Parameters:

x (np.ndarray, optional) – Final solution vector.
n_it (int, optional, default is 0) – Number of iterations performed.
F (np.ndarray, optional) – Function evaluation at the solution.
F_norm (float, optional) – Euclidean norm of F at the solution.
J (np.ndarray, optional) – Jacobian matrix at the solution.
code (TerminationCode, optional) – Termination code. Valid codes and their values are given by the TerminationCode enum.
text (str, optional) – Human-readable termination description.
success (bool, optional, default is False) – True if the solver converged.
store_iterates (bool, optional, default is False) – If True, iteration history is stored in an iterates attribute. The Iterates instance contains lists of all iterates x, function evaluations F, and step lengths relative to the full Newton step lmda.

x: ndarray | None = None¶

n_it: int = 0¶

F: ndarray | None = None¶

F_norm: float | None = None¶

J: ndarray | None = None¶

code: TerminationCode | None = None¶

text: str | None = None¶

success: bool = False¶

store_iterates: bool = False¶

terminate(converged: bool, minimum_lmda: bool, persistent_bound_violation: bool, lmda_bounds: tuple[float, float], n_it: int, max_iterations: int, violated_constraints: list[int], singular_jacobian: bool, condition_number: float) → None[source]¶: Set the termination status of the solver.

class burnman.optimize.nonlinear_solvers.DampedNewtonSolver(F, J, guess, tol: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] = 1e-06, F_tol: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] = 1e-08, max_iterations: int = 100, lambda_bounds=<function DampedNewtonSolver.<lambda>>, linear_constraints=(0.0, array([-1.])), store_iterates: bool = False, regularization: float = np.float64(2.220446049250313e-16), cond_lu_thresh: float = 1000000000000.0, constraint_thresh: float = np.float64(4.440892098500626e-16))[source]¶

Bases: object

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).

This implementation follows the algorithm 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 a tuple argument named linear_constraints. The constraints are satisfied if A*x + b <= 0. Note that the sign of b is opposite to that used in scipy.optimize.LinearConstraint. If any constraints are violated by the current candidate solution using step size lambda, lambda is reduced to satisfy all constraints.

If the current iterate lies on one or more constraints and the Newton step violates one or more of those constraints, 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 requires meeting all of the following criteria:

all(np.abs(dx (simplified Newton step)) < tol)

all(np.abs(dx (full Newton step)) < np.sqrt(10*tol))
(to avoid pathological cases)

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:

lambda is reduced to its minimum value (indicating a very nonlinear problem)

successive iterations have descent vectors which violate the constraints

the maximum number of iterations is reached

solve() → Solution[source]¶

Execute the damped Newton solver to find a root of F(x) = 0, optionally subject to linear inequality constraints A.x + b <= 0.

If the solver fails to converge, it terminates with an informative status code and description.

Returns:: Solution object. See Solution class for details of its attributes.
Return type:: Solution instance

burnman.optimize.nonlinear_solvers.damped_newton_solve(F, J, guess, tol: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] = 1e-06, F_tol: float | ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] = 1e-08, max_iterations: int = 100, lambda_bounds=<function <lambda>>, linear_constraints=(0.0, array([-1.])), store_iterates: bool = False) → Solution[source]¶

Helper function to run the DampedNewtonSolver.

Parameters:

F (callable[[np.ndarray], np.ndarray]) – Function that evaluates the system of nonlinear equations F(x) = 0 at a given x.
J (callable[[np.ndarray], np.ndarray]) – Function that evaluates the Jacobian matrix J(x) of F(x) at a given x.
guess (np.ndarray) – Initial guess for the solution vector x.
tol (float, optional) – Convergence tolerance for Newton iterations, defaults to 1.0e-6.
max_iterations (int, optional) – Maximum number of Newton iterations to perform, defaults to 100.
lambda_bounds (callable[[np.ndarray, np.ndarray], tuple[float, float]], optional) – Function returning (min_lambda, max_lambda) for the damping factor given the current Newton step dx and point x, defaults to lambda dx, x: (1.0e-8, 1.0).
linear_constraints (tuple[np.ndarray, np.ndarray], optional) – Tuple (A, b) representing linear inequality constraints A.x + b <= 0, defaults to (0.0, np.array([-1.0])).
store_iterates (bool, optional) – If True, stores all intermediate iterates, function evaluations, and lambda values in the solution object, defaults to False.

Returns:

Solution object. See Solution class for details of its attributes.

Return type:

Solution instance