4.4. mosaic_tools.ply Documentation

4.4.1. ArrayTools class

class twinkle.utils.mosaic_tools.ArrayTools

Bases: object

A few functions to work with intersection and union of numpy arrays

dict2list(Dict, keys, extra='')

to turn dict into a list. Dictionary values which will be converted to numpy array in order found in keys. dict: Dictionary to be converted to list keys: List or array with string values pointing to keys in dict Length of keys need not equal length of dict, but len(keys)<= len(dict)

dictExtract(Dict, keys, keySuffix='', newkeySuffix='')

To extract certain values given the input keys from the input dictionary and returns a new dictionary dict:Dictionary to be sampled from for new dictionary keys: keys associated with the input dictionary. keySuffix: suffix to be added to each key to access it newkeySuffix: if a new key suffix is to be added. Otherwise regular keys will be used

intersect_arrays(arrays)

Find intersection of all the arrays in “arrays” Returns the sorted, unique values that are in both of the input arrays

union_arrays(arrays)

Find UNION of all the arrays in “arrays” Returns the unique, sorted arrays of values that are in either of the two input arrays

4.4.2. FittingTools class

class twinkle.utils.mosaic_tools.FittingTools

Bases: object

A collection of utility functions for data fitting tasks.

Gauss2d_circle(x, p0=None)

Fits a circular 2D Gaussian centered at (0, 0).

Parameters:

• x (numpy.ndarray) – The 2D data points (x, y).

• p0 (list, optional) – Parameters for the Gaussian model, including amplitude and sigma.

Returns:

The fitted Gaussian values.

Return type:

numpy.ndarray

Gauss_fit(x, p0=None)

Fits a Gaussian function to the data.

Parameters:

• x (numpy.ndarray) – The input data.

• p0 (list, optional) – The parameters for the Gaussian (amplitude, mean, sigma).

Returns:

The fitted Gaussian values.

Return type:

numpy.ndarray

deviates_from_model(p0=None, fjac=True, x=None, y=None, err=None, func=None, logx=None, logy=None, loglog=None, **kwargs)

Returns the deviations (residuals) calculated from an input model function. This method is intended to be used by the Levenberg-Marquardt technique in the “mpfit.py” module, originally written in IDL by Mark Rivers and Sergey Koposov.

Parameters:

• p0 (list) – Initial parameters for the model to be fit.

• fjac (bool) – Flag for partial derivative calculation. Refer to MPFIT.py for details.

• x (numpy.ndarray) – Input x data (independent variable).

• y (numpy.ndarray) – Input y data (dependent variable).

• err (numpy.ndarray, optional) – Uncertainty in the y data.

• func (callable) – The model function that calculates y from x and parameters.

• logx (bool, optional) – Flag for applying a logarithmic transformation to x.

• logy (bool, optional) – Flag for applying a logarithmic transformation to y.

• loglog (bool, optional) – Flag for applying a logarithmic transformation to both x and y.

• **kwargs – Additional parameters to pass to the model function.

Returns:

A list containing:

• status (int): Status of the fit, used by the MPFIT.py module.

• residuals (numpy.ndarray): The residuals, either weighted or unweighted depending on the error input.

Return type:

list

erf_fit(x, p0=None)

Fits an error function to the data.

Parameters:

• x (numpy.ndarray) – The input data.

• p0 (list, optional) – Parameters for the error function (amplitude, mean, sigma).

Returns:

The fitted error function values.

Return type:

numpy.ndarray

exp_fit(x, A, b, c)

Fits an exponential function to the data.

Parameters:

• x (numpy.ndarray) – The input data.

• A (float) – The amplitude.

• b (float) – The rate constant.

• c (float) – The offset.

Returns:

The fitted exponential values.

Return type:

numpy.ndarray

get_InitParams(lenp)

Generates an initial set of starting parameters for a fit, depending on

the number of free parameters in the model.

Parameters:

lenp (int) – The number of free parameters in the model.

Returns:

A list of initial parameters, each set to 0.1.

Return type:

list

poly_nfit(x, p)

Determines the sampled values for a polynomial function whose order is

determined by the length of the input parameter array.

Parameters:

• x (numpy.ndarray) – Vector of sampling points.

• p (numpy.ndarray) – Vector of parameters for the polynomial (a0, a1, …, an).

Returns:

Sampled values of the polynomial function.

Return type:

numpy.ndarray

print_poly(polyn, xlabel='x', ylabel='y(x)', numformat='%.3f', coeff=None)

Returns a formatted string for a polynomial function of order N, including x and y labels. Optionally, the coefficients can be used to populate the polynomial’s parameters.

Parameters:

• polyn (int) – The order of the polynomial.

• xlabel (str, optional) – Label for the x-axis. Defaults to ‘x’.

• ylabel (str, optional) – Label for the y-axis. Defaults to ‘y(x)’.

• numformat (str, optional) – String format for displaying the coefficients. Defaults to ‘%.3f’.

• coeff (list, optional) – Coefficients for the polynomial.

Returns:

A string representation of the polynomial.

Return type:

str

resample_model(lam, flx, start, end, maxdelta=100.0, pband=None)

This will resample the input model spectrum to the specified resolution between the wavelengths input. If a filter is given, information from the filter will be used to supplement the resampling.

Parameters:

• lam (arr) – array of the wavelength (x-vals) for the spectrum to be resampled

• flx (arr) – array of the flux (y-vals) for the spectrum to be resampled

• start (float) – wavelength bounds for which the resampling should be conducted

• end (float) – wavelength bounds for which the resampling should be conducted

• maxdelta (float) – maximum difference between wavelegnths tolerated

• pband (obj) – passband object

Returns:

resampled spectra

resample_spectrum(dataSet1, dataSet2, resample2lowR=False)

Resamples two input spectra to the same wavelength scale over their common range.

Parameters:

• dataSet1 (tuple) – Each containing (wavelength, flux) arrays.

• dataSet2 (tuple) – Each containing (wavelength, flux) arrays.

• resample2lowR (bool, optional) – If True, resamples the high-resolution spectrum to the lower one.

Returns:

A tuple of resampled spectra.

Return type:

tuple

twoD_Gaussian(x, y, amplitude, xo, yo, sigma_x, sigma_y, theta, offset)

2D Gaussian function for modeling data, adapted from StackOverflow users ali_m and Kokomoking.

Parameters:

• x (numpy.ndarray) – The x-coordinate data.

• y (numpy.ndarray) – The y-coordinate data.

• amplitude (float) – The amplitude of the Gaussian.

• xo (float) – The x-coordinate of the center.

• yo (float) – The y-coordinate of the center.

• sigma_x (float) – The standard deviation along the x-axis.

• sigma_y (float) – The standard deviation along the y-axis.

• theta (float) – The rotation angle of the Gaussian.

• offset (float) – The offset value.

Returns:

The 2D Gaussian values at the specified coordinates.

Return type:

numpy.ndarray

4.4.3. PlottingTools class

class twinkle.utils.mosaic_tools.PlottingTools

Bases: object

coloraxes(ax, color)

plot_setup(axis, gridon=False, minortickson=True, ticklabel_fontsize=20, majortick_width=2.5, minortick_width=1.2, majortick_size=8, minortick_size=5, axes_linewidth=1.5, ytick_direction='in', xtick_direction='in', yaxis_right=False, ylog=False, xlog=False, bold=False, adjust_plot=True)

Changes the boring default matplotlib plotting canvas so that it looks nice and neat with thicker borders and larger tick marks as well as larger fontsizes for the axis labels. Options exist to include or exclude the plot grid and minortick mark labels – set up as boolean variables

simpleaxis1(ax)

This little tool erases the right and top axis lines

simpleaxis2(ax)

This little tool erases the botom and left axis lines

triple_axes_dist(xlog=False, xlabel='x', ylabel='y')

Sets up plots with 3 axes – one in center, one on right and one above center

plot. The purpose is to have a scatter plot or w/e in the center, and two distribution plots of the x and y parameters in the side plots.

Parameters:

• axScatter – axis object for center scatter plot

• ylog – booleans to indicate whether x and y axes of center plot are in log scale (base 10)

• xlog – booleans to indicate whether x and y axes of center plot are in log scale (base 10)

• xlabel – labels for x and y axes labels

• ylabel – labels for x and y axes labels

• Return

zeroaxes(ax)

4.4.4. ReadWrite_Tools class

class twinkle.utils.mosaic_tools.ReadWrite_Tools

Bases: object

Read and Write tools once stuff is read

create_datadict(hnames, data)

Create a dictionary from a 2D array of data and corresponding header names.

Parameters:

• hnames (list) – List of header names to use as keys.

• data (array) – 2D array where each column corresponds to a header.

Returns:

Dictionary with header names as keys and data columns as values.

Return type:

dict

create_header(list0, more=None, nowrite=None, delimiter='\t ')

Create a header string from a list of strings, excluding specified items.

Parameters:

• list0 (list) – List of strings to include in the header.

• more (list, optional) – Additional strings to append to the header.

• nowrite (list, optional) – Items to exclude from the header.

• delimiter (str, optional) – Delimiter used to separate header items. Defaults to ‘t’.

Returns:

Formatted header string with specified items and delimiter.

Return type:

str

sort_duplicates(file, dupcol='object_u', duplicates=None)

Sort duplicate values in a data file and save a new version with duplicates at the top.

Parameters:

• file (str) – Path to the input file.

• dupcol (str, optional) – Column name used to identify duplicates. Defaults to ‘object_u’.

• duplicates (list, optional) – Known duplicates. If None, duplicates are identified automatically.

Returns:

Filename of the new sorted file.

Return type:

str

4.4.5. mpfit class

class twinkle.utils.mosaic_tools.mpfit(fcn, xall=None, functkw={}, parinfo=None, ftol=1e-10, xtol=1e-10, gtol=1e-10, damp=0.0, maxiter=2000, factor=100.0, nprint=1, iterfunct='default', iterkw={}, nocovar=0, rescale=0, autoderivative=1, quiet=0, diag=None, epsfcn=None, debug=0)

Bases: object

Perform Levenberg-Marquardt least-squares minimization, based on MINPACK-1.

AUTHORS The original version of this software, called LMFIT, was written in FORTRAN as part of the MINPACK-1 package by XXX.

Craig Markwardt converted the FORTRAN code to IDL. The information for the IDL version is: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html

Mark Rivers created this Python version from Craig’s IDL version. Mark Rivers, University of Chicago Building 434A, Argonne National Laboratory 9700 South Cass Avenue, Argonne, IL 60439 rivers@cars.uchicago.edu Updated versions can be found at http://cars.uchicago.edu/software

Sergey Koposov converted the Mark’s Python version from Numeric to numpy Sergey Koposov, University of Cambridge, Institute of Astronomy, Madingley road, CB3 0HA, Cambridge, UK koposov@ast.cam.ac.uk Updated versions can be found at http://code.google.com/p/astrolibpy/source/browse/trunk/

DESCRIPTION

MPFIT uses the Levenberg-Marquardt technique to solve the least-squares problem. In its typical use, MPFIT will be used to fit a user-supplied function (the “model”) to user-supplied data points (the “data”) by adjusting a set of parameters. MPFIT is based upon MINPACK-1 (LMDIF.F) by More’ and collaborators.

For example, a researcher may think that a set of observed data points is best modelled with a Gaussian curve. A Gaussian curve is parameterized by its mean, standard deviation and normalization. MPFIT will, within certain constraints, find the set of parameters which best fits the data. The fit is “best” in the least-squares sense; that is, the sum of the weighted squared differences between the model and data is minimized.

The Levenberg-Marquardt technique is a particular strategy for iteratively searching for the best fit. This particular implementation is drawn from MINPACK-1 (see NETLIB), and is much faster and more accurate than the version provided in the Scientific Python package in Scientific.Functions.LeastSquares. This version allows upper and lower bounding constraints to be placed on each parameter, or the parameter can be held fixed.

The user-supplied Python function should return an array of weighted deviations between model and data. In a typical scientific problem the residuals should be weighted so that each deviate has a gaussian sigma of 1.0. If X represents values of the independent variable, Y represents a measurement for each value of X, and ERR represents the error in the measurements, then the deviates could be calculated as follows:

DEVIATES = (Y - F(X)) / ERR

where F is the analytical function representing the model. You are recommended to use the convenience functions MPFITFUN and MPFITEXPR, which are driver functions that calculate the deviates for you. If ERR are the 1-sigma uncertainties in Y, then

TOTAL( DEVIATES^2 )

will be the total chi-squared value. MPFIT will minimize the chi-square value. The values of X, Y and ERR are passed through MPFIT to the user-supplied function via the FUNCTKW keyword.

Simple constraints can be placed on parameter values by using the PARINFO keyword to MPFIT. See below for a description of this keyword.

MPFIT does not perform more general optimization tasks. See TNMIN instead. MPFIT is customized, based on MINPACK-1, to the least-squares minimization problem.

USER FUNCTION

The user must define a function which returns the appropriate values as specified above. The function should return the weighted deviations between the model and the data. It should also return a status flag and an optional partial derivative array. For applications which use finite-difference derivatives – the default – the user function should be declared in the following way:

def myfunct(p, fjac=None, x=None, y=None, err=None) # Parameter values are passed in “p” # If fjac==None then partial derivatives should not be # computed. It will always be None if MPFIT is called with default # flag. model = F(x, p) # Non-negative status value means MPFIT should continue, negative means # stop the calculation. status = 0 return([status, (y-model)/err]

See below for applications with analytical derivatives.

The keyword parameters X, Y, and ERR in the example above are suggestive but not required. Any parameters can be passed to MYFUNCT by using the functkw keyword to MPFIT. Use MPFITFUN and MPFITEXPR if you need ideas on how to do that. The function must accept a parameter list, P.

In general there are no restrictions on the number of dimensions in X, Y or ERR. However the deviates must be returned in a one-dimensional Numeric array of type Float.

User functions may also indicate a fatal error condition using the status return described above. If status is set to a number between -15 and -1 then MPFIT will stop the calculation and return to the caller.

ANALYTIC DERIVATIVES

In the search for the best-fit solution, MPFIT by default calculates derivatives numerically via a finite difference approximation. The user-supplied function need not calculate the derivatives explicitly. However, if you desire to compute them analytically, then the AUTODERIVATIVE=0 keyword must be passed to MPFIT. As a practical matter, it is often sufficient and even faster to allow MPFIT to calculate the derivatives numerically, and so AUTODERIVATIVE=0 is not necessary.

If AUTODERIVATIVE=0 is used then the user function must check the parameter FJAC, and if FJAC!=None then return the partial derivative array in the return list.

where FGRAD(x, p, i) is a user function which must compute the derivative of the model with respect to parameter P[i] at X. When finite differencing is used for computing derivatives (ie, when AUTODERIVATIVE=1), or when MPFIT needs only the errors but not the derivatives the parameter FJAC=None.

Derivatives should be returned in the PDERIV array. PDERIV should be an m x n array, where m is the number of data points and n is the number of parameters. dp[i,j] is the derivative at the ith point with respect to the jth parameter.

The derivatives with respect to fixed parameters are ignored; zero is an appropriate value to insert for those derivatives. Upon input to the user function, FJAC is set to a vector with the same length as P, with a value of 1 for a parameter which is free, and a value of zero for a parameter which is fixed (and hence no derivative needs to be calculated).

If the data is higher than one dimensional, then the last dimension should be the parameter dimension. Example: fitting a 50x50 image, “dp” should be 50x50xNPAR.

CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD

The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another.

These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT.

PARINFO should be a list of dictionaries, one list entry for each parameter. Each parameter is associated with one element of the array, in numerical order. The dictionary can have the following keys (none are required, keys are case insensitive):

‘value’ - the starting parameter value (but see the START_PARAMS

parameter for more information).

‘fixed’ - a boolean value, whether the parameter is to be held

fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation.

‘limited’ - a two-element boolean array. If the first/second

element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together.

‘limits’ - a two-element float array. Gives the

parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together.

‘parname’ - a string, giving the name of the parameter. The

fitting code of MPFIT does not use this tag in any way. However, the default iterfunct will print the parameter name if available.

‘step’ - the step size to be used in calculating the numerical

derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0.

‘mpside’ - the sidedness of the finite difference when computing

numerical derivatives. This field can take four values:

0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h

-1 - one-sided derivative (f(x) - f(x-h))/h

2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)

Where H is the STEP parameter described above. The “automatic” one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0.

‘mpmaxstep’ - the maximum change to be made in the parameter

value. During the fitting process, the parameter will never be changed by more than this value in one iteration.

A value of 0 indicates no maximum. Default: 0.

‘tied’ - a string expression which “ties” the parameter to other

free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo(2).tied = ‘2 * p(1)’. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can’t be used in expressions. ]

‘mpprint’ - if set to 1, then the default iterfunct will print the

parameter value. If set to 0, the parameter value will not be printed. This tag can be used to selectively print only a few parameter values out of many. Default: 1 (all parameters printed)

Future modifications to the PARINFO structure, if any, will involve adding dictionary tags beginning with the two letters “MP”. Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored.

PARINFO Example: parinfo = [{‘value’:0., ‘fixed’:0, ‘limited’:[0,0], ‘limits’:[0.,0.]}

for i in range(5)]

parinfo[0][‘fixed’] = 1 parinfo[4][‘limited’][0] = 1 parinfo[4][‘limits’][0] = 50. values = [5.7, 2.2, 500., 1.5, 2000.] for i in range(5): parinfo[i][‘value’]=values[i]

A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50.

EXAMPLE

import mpfit import np.oldnumeric as Numeric x = arange(100, float) p0 = [5.7, 2.2, 500., 1.5, 2000.] y = ( p[0] + p[1]*[x] + p[2]*[x**2] + p[3]*sqrt(x) + p[4]*log(x)) fa = {‘x’:x, ‘y’:y, ‘err’:err} m = mpfit(‘myfunct’, p0, functkw=fa) print ‘status = ‘, m.status if (m.status <= 0): print ‘error message = ‘, m.errmsg print ‘parameters = ‘, m.params

Minimizes sum of squares of MYFUNCT. MYFUNCT is called with the X, Y, and ERR keyword parameters that are given by FUNCTKW. The results can be obtained from the returned object m.

THEORY OF OPERATION

There are many specific strategies for function minimization. One very popular technique is to use function gradient information to realize the local structure of the function. Near a local minimum the function value can be taylor expanded about x0 as follows:

f(x) = f(x0) + f’(x0) . (x-x0) + (1/2) (x-x0) . f’’(x0) . (x-x0) Order 0th 1st 2nd

Here f’(x) is the gradient vector of f at x, and f’’(x) is the Hessian matrix of second derivatives of f at x. The vector x is the set of function parameters, not the measured data vector. One can find the minimum of f, f(xm) using Newton’s method, and arrives at the following linear equation:

f’’(x0) . (xm-x0) = - f’(x0) (2)

If an inverse can be found for f’’(x0) then one can solve for (xm-x0), the step vector from the current position x0 to the new projected minimum. Here the problem has been linearized (ie, the gradient information is known to first order). f’’(x0) is symmetric n x n matrix, and should be positive definite.

The Levenberg - Marquardt technique is a variation on this theme. It adds an additional diagonal term to the equation which may aid the convergence properties:

(f’’(x0) + nu I) . (xm-x0) = -f’(x0) (2a)

where I is the identity matrix. When nu is large, the overall matrix is diagonally dominant, and the iterations follow steepest descent. When nu is small, the iterations are quadratically convergent.

In principle, if f’’(x0) and f’(x0) are known then xm-x0 can be determined. However the Hessian matrix is often difficult or impossible to compute. The gradient f’(x0) may be easier to compute, if even by finite difference techniques. So-called quasi-Newton techniques attempt to successively estimate f’’(x0) by building up gradient information as the iterations proceed.

In the least squares problem there are further simplifications which assist in solving eqn (2). The function to be minimized is a sum of squares:

f = Sum(hi^2) (3)

where hi is the ith residual out of m residuals as described above. This can be substituted back into eqn (2) after computing the derivatives:

f’ = 2 Sum(hi hi’) f’’ = 2 Sum(hi’ hj’) + 2 Sum(hi hi’’) (4)

If one assumes that the parameters are already close enough to a minimum, then one typically finds that the second term in f’’ is negligible [or, in any case, is too difficult to compute]. Thus, equation (2) can be solved, at least approximately, using only gradient information.

In matrix notation, the combination of eqns (2) and (4) becomes:

hT’ . h’ . dx = - hT’ . h (5)

Where h is the residual vector (length m), hT is its transpose, h’ is the Jacobian matrix (dimensions n x m), and dx is (xm-x0). The user function supplies the residual vector h, and in some cases h’ when it is not found by finite differences (see MPFIT_FDJAC2, which finds h and hT’). Even if dx is not the best absolute step to take, it does provide a good estimate of the best direction, so often a line minimization will occur along the dx vector direction.

The method of solution employed by MINPACK is to form the Q . R factorization of h’, where Q is an orthogonal matrix such that QT . Q = I, and R is upper right triangular. Using h’ = Q . R and the ortogonality of Q, eqn (5) becomes

(RT . QT) . (Q . R) . dx = - (RT . QT) . h RT . R . dx = - RT . QT . h (6) R . dx = - QT . h

where the last statement follows because R is upper triangular. Here, R, QT and h are known so this is a matter of solving for dx. The routine MPFIT_QRFAC provides the QR factorization of h, with pivoting, and MPFIT_QRSOLV provides the solution for dx.

REFERENCES

MINPACK-1, Jorge More’, available from netlib (www.netlib.org). “Optimization Software Guide,” Jorge More’ and Stephen Wright, SIAM, Frontiers in Applied Mathematics, Number 14. More’, Jorge J., “The Levenberg-Marquardt Algorithm: Implementation and Theory,” in Numerical Analysis, ed. Watson, G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.

MODIFICATION HISTORY

Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM Copyright (C) 1997-2002, Craig Markwardt This software is provided as is without any warranty whatsoever. Permission to use, copy, modify, and distribute modified or unmodified copies is granted, provided this copyright and disclaimer are included unchanged.

Translated from MPFIT (Craig Markwardt’s IDL package) to Python, August, 2002. Mark Rivers Converted from Numeric to numpy (Sergey Koposov, July 2008)

blas_enorm32 = <fortran function dnrm2>

blas_enorm64 = <fortran function dnrm2>

calc_covar(rr, ipvt=None, tol=1e-14)

Calculate the covariance matrix based on the QR factorization.

Parameters:

• rr (array) – The upper triangular matrix from QR factorization.

• ipvt (array, optional) – Pivot indices.

• tol (float, optional) – Tolerance for determining if a parameter is significant.

Returns:

The covariance matrix or -1 in case of an error.

Return type:

array

call(fcn, x, functkw, fjac=None)

Call the user-defined function with the provided parameters.

Parameters:

• fcn (callable) – The user’s function to be called.

• x (array) – The current parameter estimates to be passed to fcn.

• functkw (dict) – Additional keyword arguments to be passed to fcn.

• fjac (optional) – Jacobian matrix or related parameter information.

Returns:

A list containing:

• status (int): Status of the function call.

• f (array): Function values calculated at x.

Return type:

list

defiter(fcn, x, iter, fnorm=None, functkw=None, quiet=0, iterstop=None, parinfo=None, format=None, pformat='%.10g', dof=1)

Print the current iteration information during optimization.

Parameters:

• fcn (callable) – The function to be minimized.

• x (array) – Current parameter estimates.

• iter (int) – The current iteration number.

• fnorm (float, optional) – Current function norm (chi-squared). If not provided, it will be computed.

• functkw (dict, optional) – Additional keyword arguments to be passed to fcn.

• quiet (int, optional) – If set to a non-zero value, suppress output. Default is 0 (output shown).

• iterstop (optional) – A hypothetical parameter for iteration stopping condition.

• parinfo (list of dicts, optional) – Information for each parameter such as names.

• format (optional) – Custom format for printing.

• pformat (str, optional) – Format for printing parameters. Default is ‘%.10g’.

• dof (int, optional) – Degrees of freedom. Default is 1.

Returns:

Always returns 0 after printing.

Return type:

int

enorm(vec)

Compute the Euclidean norm of a vector.

Parameters:

vec (array) – The input vector.

Returns:

The Euclidean norm of the input vector.

Return type:

float

fdjac2(fcn, x, fvec, step=None, ulimited=None, ulimit=None, dside=None, epsfcn=None, autoderivative=1, functkw=None, xall=None, ifree=None, dstep=None)

Compute the Jacobian matrix using finite differences.

Parameters:

• fcn (callable) – The function to compute the Jacobian for.

• x (array) – The current parameter estimates.

• fvec (array) – Function values at x.

• step (array, optional) – Steps for finite differencing. Default is None.

• ulimited (array, optional) – Flags indicating whether parameters are unlimited.

• ulimit (array, optional) – Upper limits for parameters.

• dside (array, optional) – Side of finite difference to use.

• epsfcn (float, optional) – Parameter for finite differencing.

• autoderivative (int, optional) – Set to 0 to use analytical derivatives instead.

• functkw (dict, optional) – Additional keyword arguments for fcn.

• xall (array, optional) – Array of all parameters.

• ifree (array, optional) – Indices of free parameters.

• dstep (array, optional) – Absolute or relative step size.

Returns:

Jacobian matrix evaluated at x.

Return type:

array

lmpar(r, ipvt, diag, qtb, delta, x, sdiag, par=None)

Perform damped least squares minimization.

Parameters:

• r (array) – Upper triangular matrix from QR factorization.

• ipvt (array) – Pivot indices.

• diag (array) – Diagonal scaling factors.

• qtb (array) – The product of (Q^T) and the right-hand side vector b.

• delta (float) – Convergence threshold.

• x (array) – Initial guess for the solution.

• sdiag (array) – Storage for diagonal elements from the QR factorization.

• par (float, optional) – Initial value for the parameter.

Returns:

A list containing:

• r (array): Modified upper triangular matrix.

• par (float): Updated parameter value.

• x (array): Solution vector.

• sdiag (array): Updated diagonal elements.

Return type:

list

parinfo(parinfo=None, key='a', default=None, n=0)

Extract values from a parameter information structure.

Parameters:

• parinfo (list of dicts, optional) – List containing parameter information.

• key (str, optional) – The key to look for within parinfo. Default is ‘a’.

• default (optional) – Default value to return if the key is not found.

• n (int, optional) – The number of parameters to return. Defaults to 0.

Returns:

A list of values extracted from parinfo, or default values if not found.

Return type:

list

qrfac(a, pivot=0)

Perform QR factorization of a matrix.

Parameters:

• a (array) – Input matrix to be factorized.

• pivot (int, optional) – Whether to perform pivoting (0 for no, 1 for yes).

Returns:

A list containing:

• a (array): The upper triangular matrix after factorization.

• ipvt (array): The pivot indices.

• rdiag (array): The diagonal elements of the R matrix.

• acnorm (array): The norms of the columns of A.

Return type:

list

qrsolv(r, ipvt, diag, qtb, sdiag)

Solve a linear system using QR factorization.

Parameters:

• r (array) – Upper triangular matrix from QR factorization.

• ipvt (array) – Pivot indices.

• diag (array) – Diagonal scaling factors.

• qtb (array) – The product of (Q^T) and the right-hand side vector b.

• sdiag (array) – Storage for diagonal elements.

Returns:

A tuple containing:

• r (array): Modified upper triangular matrix.

• x (array): Solution vector.

• sdiag (array): Updated diagonal elements.

Return type:

tuple

tie(p, ptied=None)

Apply constraints by tying parameters to each other.

Parameters:

• p (array) – Current parameters.

• ptied (list of str, optional) – List of parameter names to tie.

Returns:

The modified parameters.

Return type:

array