pyextal.jaxMultipole.fit ======================== .. py:module:: pyextal.jaxMultipole.fit .. autoapi-nested-parse:: Weighted least-squares fit of a multipole model to structure factors. The objective is the usual .. math:: \chi^2 = \sum_i \left(\frac{y_i^{\rm obs} - y_i^{\rm calc}(\mathbf x)}{\sigma_i}\right)^2 + \sum_r \left(\frac{g_r(\mathbf p) - t_r}{\sigma_r}\right)^2 over the observations (real parts, and imaginary parts where an observation carries one) and any soft restraints. It is minimised with ``scipy.optimize.least_squares`` (trust-region reflective, box bounds), with an exact Jacobian from ``jax.jacfwd`` -- in float64, inside ``jax.enable_x64``, so the rest of the process (the float32 Bloch engine) is unaffected. What is reported ---------------- * ``esd`` -- from ``(J^T J)^-1``, i.e. taking the sigmas at face value, and ``esd_scaled`` -- the same times the goodness of fit, the conventional crystallographic choice when GoF > 1; * the correlation matrix of the free parameters and the pairs above 0.9; * ``R = sum |dy| / sum |y|``, ``wR = sqrt(sum w dy^2 / sum w y^2)``, GoF ``= sqrt(chi^2 / (n - m))``; * the singular values and condition number of the weighted Jacobian. :func:`design_check` gives the last two *before* fitting, at the start point: which directions in parameter space the data do not constrain, and which parameters they load on. With ~25 reflections and ADPs refined it is the first thing to look at. :func:`hamilton_test` compares nested models. Device: the fit runs on the CPU by default (``device='cpu'``). With some 30 reflections it is a few hundred small float64 kernels, which a CPU does in ~15 s for the STO ladder while a consumer GPU -- float64 at a small fraction of its float32 rate, plus launch overhead -- takes ~50 s, to identical results. ``device='default'`` uses whatever JAX picks; the GPU matters for refining against intensities (Option 2), not here. Classes ------- .. autoapisummary:: pyextal.jaxMultipole.fit.Restraint pyextal.jaxMultipole.fit.FitResult pyextal.jaxMultipole.fit.Fitter Functions --------- .. autoapisummary:: pyextal.jaxMultipole.fit.residual_fn pyextal.jaxMultipole.fit.fit pyextal.jaxMultipole.fit.design_check pyextal.jaxMultipole.fit.hamilton_test Module Contents --------------- .. py:class:: Restraint A soft restraint on one physical parameter: ``(p[name] - target) / sigma``. .. py:attribute:: name :type: str .. py:attribute:: target :type: float .. py:attribute:: sigma :type: float .. py:class:: FitResult Outcome of :func:`fit`. Attributes are NumPy arrays or floats; ``names`` label ``x``, ``esd`` and the rows/columns of ``cov``/``corr``; ``p`` and ``p_esd`` follow the model layout. .. py:attribute:: names :type: list .. py:attribute:: x :type: numpy.ndarray .. py:attribute:: p :type: numpy.ndarray .. py:attribute:: p_esd :type: numpy.ndarray .. py:attribute:: cov :type: numpy.ndarray .. py:attribute:: esd :type: numpy.ndarray .. py:attribute:: esd_scaled :type: numpy.ndarray .. py:attribute:: corr :type: numpy.ndarray .. py:attribute:: chi2 :type: float .. py:attribute:: n_obs :type: int .. py:attribute:: n_par :type: int .. py:attribute:: gof :type: float .. py:attribute:: R :type: float .. py:attribute:: wR :type: float .. py:attribute:: calc :type: numpy.ndarray .. py:attribute:: obs :type: object .. py:attribute:: singular_values :type: numpy.ndarray :value: None .. py:attribute:: condition :type: float :value: None .. py:attribute:: success :type: bool :value: True .. py:attribute:: message :type: str :value: '' .. py:attribute:: pmap :type: object :value: None .. py:method:: value(name) .. py:method:: p_esd_of(name) .. py:method:: has(name) .. py:method:: strong_correlations(threshold=0.9) .. py:function:: residual_fn(model, obs, pmap, restraints=()) ``x -> weighted residual vector`` (real), plus ``x -> calc``. Real parts always count; imaginary parts only where ``obs.use_imag``. .. py:class:: Fitter(model, obs, restraints=(), device='cpu') A model and a data set with their residual and Jacobian compiled once. The residuals are compiled as a function of the *full* physical vector p, so any `ParameterMap` -- any choice of frozen and free parameters -- reuses them: ``r(x) = r_p(p0 + M x)``, ``J(x) = J_p(p) M``. Refinements that change what is refined many times over (ladders, iterative block schemes, chi^2 profiles) then cost one compilation instead of one per fit. `fit` keeps a small cache of these, so ordinary calls benefit too. .. py:attribute:: restraints :value: () .. py:method:: fit(pmap, x0=None, max_nfev=200, verbose=0) Fits the free parameters of ``pmap``; see :func:`fit`. .. py:function:: fit(model, obs, pmap, restraints=(), x0=None, max_nfev=200, verbose=0, device='cpu') Fits the free parameters of ``pmap`` to ``obs``. :param model: The model. :type model: MultipoleModel :param obs: Data (`jaxMultipole.data`). :type obs: Observations :param pmap: Free parameters (`jaxMultipole.params`). :type pmap: ParameterMap :param restraints: Soft restraints. :type restraints: list[Restraint] :param x0: Start; ``pmap.x0`` by default. :type x0: array, optional :param max_nfev: Function-evaluation budget. :type max_nfev: int :param verbose: ``least_squares`` verbosity. :type verbose: int :param device: ``'cpu'`` (default), ``'gpu'`` or ``'default'``. :type device: str :returns: FitResult .. note:: The compiled residuals are cached per (model, obs, restraints, device) object, so repeated fits of the same data with different free parameters compile once (`Fitter`). .. py:function:: design_check(model, obs, pmap, restraints=(), x=None, top=4, device='cpu') How well the data constrain each direction in parameter space. :returns: ``singular_values``, ``condition``, predicted ``esd`` and ``corr`` at ``x`` (the start by default), and ``weak`` -- for each singular value below 1e-3 of the largest, the parameters it loads on. :rtype: dict .. py:function:: hamilton_test(fit_restricted, fit_full, alpha=0.05) Hamilton's R-factor ratio test of a nested pair of fits. :param fit_restricted: The model with fewer parameters. :type fit_restricted: FitResult :param fit_full: The model with ``b`` more. :type fit_full: FitResult :param alpha: Significance level. :type alpha: float :returns: ``ratio = wR_restricted / wR_full``, the critical value ``sqrt(1 + b/(n-m) F_{b, n-m, alpha})`` and ``significant``. :rtype: dict