Optimization#
- class burnman.optimize.composition_fitting.DummyCompositionSolution(endmember_element_formulae, endmember_site_formulae)[source]#
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:
- 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
adiabatic_bulk_modulus()
- property K_T#
Alias for
isothermal_bulk_modulus()
- 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 adiabatic_bulk_modulus#
Returns adiabatic bulk modulus of the solution [Pa]. Aliased with self.K_S.
- property adiabatic_bulk_modulus_reuss#
Alias for
adiabatic_bulk_modulus()
- property adiabatic_compressibility#
Returns adiabatic compressibility of the solution. (or inverse adiabatic bulk modulus) [1/Pa]. Aliased with self.K_S.
- property adiabatic_compressibility_reuss#
Alias for
adiabatic_compressibility()
- property alpha#
Alias for
thermal_expansivity()
- property beta_S#
Alias for
adiabatic_compressibility()
- property beta_T#
Alias for
isothermal_compressibility()
- property bulk_sound_velocity#
Returns bulk sound speed of the solution [m/s]. Aliased with self.v_phi.
- 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)#
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:
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] for temperatures[i] and pressures[i].
- Return type:
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 isothermal_bulk_modulus#
Returns isothermal bulk modulus of the solution [Pa]. Aliased with self.K_T.
- property isothermal_bulk_modulus_reuss#
Alias for
isothermal_bulk_modulus()
- property isothermal_compressibility#
Returns isothermal compressibility of the solution. (or inverse isothermal bulk modulus) [1/Pa]. Aliased with self.K_T.
- property isothermal_compressibility_reuss#
Alias for
isothermal_compressibility()
- 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:
- 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.
- set_composition(molar_fractions)#
Set the composition for this solution. Resets cached properties.
- 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.
- 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 [5.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.
- 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:
- 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:
- 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)[source]#
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)[source]#
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)[source]#
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.
- 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)[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.fit_PTV_data(mineral, fit_params, data, data_covariances=[], delta_params=None, bounds=None, param_tolerance=1e-05, max_lm_iterations=50, verbose=True)[source]#
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)[source]#
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.
- 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)[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.linear_fitting.weighted_constrained_least_squares(A, b, Cov_b=None, equality_constraints=None, inequality_constraints=None)[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.
- 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:
- burnman.optimize.nonlinear_fitting.nonlinear_least_squares_fit(model, lm_damping=0.0, param_tolerance=1e-07, max_lm_iterations=100, 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:
- :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)[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_variancex_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:
- Returns:
matplotlib.pyplot.figure
and list ofmatplotlib.pyplot.Axes
objects.- Return type:
- 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.
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:
- 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:
- burnman.optimize.nonlinear_fitting.plot_residuals(ax, weighted_residuals, n_bins=None, flags=[])[source]#
- burnman.optimize.nonlinear_solvers.solve_constraint_lagrangian(x, jac_x, c_x, c_prime)[source]#
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:
- 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)[source]#
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