onls 0.3 (09-11-2026)

NEW FUNCTIONS
* doGauge() and demo "gauge": onls is checked in 209 comparisons: 129 reference values, of four kinds: published orthogonal-regression results (ODRPACK guide example 2.C.i, Algorithm 676, the bounded ODRPACK95 example, Pearson-York, Daeron & Vermeesch 2024, the XLSTAT Deming example), closed-form solutions (total least squares, generalized total least squares, Deming regression), ODRPACK 2.01 run to convergence (example with separate WE/WD weights; the NIST problems Chwirut2, Gauss1, ENSO and the Rat43 data as orthogonal regressions), and the NIST certified values (ordinary least squares, reproduced by onls in the limit of a negligible sigma_x), and 80 consistency checks of the functions around the fit and of the agreement of the two algorithms (replacing the testthat suite of 0.2-1). Each comparison is reported as PASS, NOTE (documented discrepancy of the reference) or DIFF; the output goes to the console, an HTML or a PDF file. Current result: 206 PASS, 3 NOTE, 0 DIFF. The reference examples 4 to 12 of ?onls moved into the demo
* anova.onls(): comparison of nested onls models by the minimized orthogonal objective S, with an extra-sum-of-squares F test for estimated variances and a likelihood-ratio chi-square test (S_a - S_b) for known variances. Before, anova() dispatched to anova.nls(), which compares the vertical residual sums of squares (deviance()), a criterion that onls does not minimize. Models fitted to different observations, with different weights, sigma_y, sigma_x or known_sigma are refused. In simulations under the null hypothesis the tests kept the 5% level (4.5-6.0% rejections, against 3.5-7.0% for the vertical test)

REMOVED FUNCTIONS
* NIST() was removed; its fits are part of the demo "gauge". The Gauss1 data in NIST() had been rounded to 4 decimals, so that their least-squares fit missed the NIST certified residual sum of squares (1315.826 instead of 1315.822); the demo uses the full-precision NIST data

UPDATED FUNCTIONS
* onls():
  - new argument `algorithm = c("auto", "structured", "dense")`. "structured" computes each Levenberg-Marquardt step of the joint problem by eliminating the p x p blocks of the foot-point corrections (Schur complement, as in ODRPACK), so that the cost per iteration grows linearly with n instead of with n^3: a fit with n = 1000 takes about 0.02 s instead of 15 s, n = 100 000 about 2 s. The damping follows MINPACK's trust-region rule, so on well-posed problems the iterates and the estimates agree with "dense" (the solver of 0.2, minpack.lm::nls.lm on the full Jacobian) to about 1e-15. The elimination is written without differences of nearly equal numbers, so that it stays accurate when the response is much more precise than the predictors (e.g. y in much larger units than x with the default unit precisions). Correlated sigma_x is solved by a Cholesky factorization vectorized over the observations. Parameter bounds are handled by an active set: parameters are held at a bound only while the descent direction points outward, so the structured algorithm reaches the box optimum, where the dense algorithm (minpack.lm's bounded LM) can stop 5-25 % above it, with non-orthogonal points. The default "auto" uses "structured" for more than 100 observations or when parameter bounds are given, and "dense" otherwise (for small n the compiled minpack.lm code has less overhead per iteration). Fits with the foot-point bounds `extend`/`window` use "dense". New example 4 in ?onls and vignette section 9.3 compare the two algorithms
  - `weights`, `subset`, `sigma_x` and `sigma_y` are evaluated in `data` first, as in nls(), so that columns can be named directly (e.g. `weights = w, subset = conc > 0.1, sigma_x = sx`); NA in a logical `subset` drops the observation
  - all forms of `sigma_x` (none, scalar, length-p, p x p covariance, n x p per observation) are stored and used as one n x p x p array of Cholesky factors of the precisions (`Lx`), instead of separate diagonal, matrix and per-observation code paths
  - the covariance matrix is inverted after equilibration (both algorithms), so that parameters on very different scales (e.g. 1e-9 and 2) no longer give "Covariance matrix singular" and NA standard errors
  - observations with zero weight no longer count towards `df_resid`
  - with all parameters `fixed`, only the foot points are optimized (used to profile one-parameter models)
  - new warnings: when the variances are treated as known (e.g. only `sigma_y` given) but `sigma_x` is missing, because the predictor errors are then taken as standard deviation 1 in the units of the predictors; and when, with `extend`/`window` foot-point bounds, the bounded solver stopped at foot points that are neither orthogonal nor at their bound. Clear errors for NA in `sigma_x` or `sigma_y`
  - the fit object was trimmed from 60 to 36 elements by removing duplicates and internal quantities. Removed: `parONLS` (use coef()), `x0`, `y0` (use x0(), y0() or `xi`, `yi`), `resid_o` (residuals_o() or `orth_dist`), `pred`, `resp` (`X_obs`, `y_obs`), `residONLS`, `residNLS`, `fittedNLS` (residuals(), fitted()), `grad`, `QR`, `x_corrections`, `Q_x`, `Sigma_x`, `is_diag_x`, `pred_names`, `param_names`, `resp_name`, `start`, `control_in`, `iters` (convInfo$finIter), `p`, `q`, `type` and `model`; `control` is now the `control` argument as given. New: `algorithm`, `Lx`, `yi` (foot-point responses), `extend` and `window`. All per-observation elements are in the original observation order. ?onls documents all elements
* logLik_o() is now the log-likelihood of the errors-in-variables model: each observation contributes 1 + p Gaussian error components (response and predictors) with the precisions of the fit; the variance scale is 1 for known variances (the value was previously computed as if the scale were estimated) and its ML estimate otherwise. "df" counts the free parameters, plus one for an estimated scale (previously always 1 + q). Observations with zero weight are excluded (it returned -Inf). Values differ from earlier versions; AIC/BIC compare onls fits to the same data with the same error variances
* predict.onls():
  - prediction intervals for a new observation of the errors-in-variables model, with variance se_fit^2 + s^2 (sigma_y^2 + g' Sigma_x g) (g: slope of the model in the predictors). Previously the prediction interval added the variance of the weighted vertical residuals, which ignores the predictor error and the known variances
  - new arguments `sigma_x` and `sigma_y` for the new observations (required if the fit used observation-specific values; `sigma_x = 0` for an exactly known predictor), and `scale` and `df`; `scale` also rescales the standard errors of the fitted values, as in predict.lm()
  - normal quantiles when the variances are known, t quantiles with df_resid otherwise
  - failed with "non-conformable arguments" for `se.fit` or intervals when parameters were `fixed`, contrary to the NEWS of 0.2
  - evaluates the model formula in its environment, so that functions and constants defined where the model was fitted are found
  - new help page ?predict.onls
* summary.onls():
  - normal quantiles when the variances are known ("z value", "Pr(>|z|)"), as confint.onls() already did, t quantiles with df_resid otherwise
  - `correlation = TRUE` failed with "non-conformable arrays" when parameters were `fixed`
* confint.onls(): new argument `method = c("profile", "wald", "boot")`; the default is now the profile interval. The bootstrap could loop forever (refits were retried without limit until `k` succeeded); it now stops after `max_tries` attempts with a warning. It resampled the data rows but not the per-observation `weights`, `sigma_x` and `sigma_y`, which stayed in the original order and were paired with the wrong rows; they are now resampled with their rows. The bootstrap no longer re-evaluates the `data` expression in the calling frame (which failed for fits made inside functions), no longer discards estimates beyond 20 standard errors, refits with the algorithm of the original fit, and prints progress only with `verbose = TRUE`. `parm` accepts names and indices
* profile.onls() is now registered as an S3 method, so profile() works (before, profile() dispatched to profile.nls and failed). It refits from the data stored in the fit, with the algorithm of the original fit, uses s2 = 1 and chi-square cut-offs for known variances, and gained plot() and confint() methods
* check_o() measures the angle between the residual vector and the model surface in coordinates whitened by the precisions. This single criterion replaces the tangent angle (unweighted) and the per-axis relative KKT residual (weighted), and fixes false "FAILED" messages for a correlated (p x p) `sigma_x`, where all points were reported as non-orthogonal, and for predictors with very small `sigma_x`. `tol_kkt` is no longer used. Values are returned in the original observation order
* x0(), y0() return values in the original observation order, like residuals_o(). With `na.action = na.exclude`, NA were inserted at the original positions into sorted values, so that values sat at the wrong rows, and the orthogonality check inside onls() compared vectors of different lengths
* plot.onls() evaluates the model formula in its environment, as predict.onls()
* print.onls() failed with "missing value where TRUE/FALSE needed" when the orthogonality check could not be computed; it now says so. The convergence line reads "Convergence tolerance (ftol)", since it shows the tolerance used, not an achieved value

onls 0.2 (21-09-2026)

MAJOR UPDATE: one-pass (ODRPACK-style) algorithm and multivariateness:

* replaced the nested two-step algorithm (inner foot-point optimization for every observation, outer optimization of the parameters, using L-BFGS-B for foot points when p > 1) by a one-pass algorithm in the manner of ODRPACK: the parameters and all foot-point corrections delta_i are estimated simultaneously by a single Levenberg-Marquardt run (minpack.lm::nls.lm) on the joint problem, minimizing the same objective as ODRPACK, sum[ Qyy_i * (y_i - f(x_i + delta_i, theta))^2 + delta_i' * Qx_i * delta_i ], with Qyy and Qx corresponding to ODRPACK's WE and WD
* the same code path now handles any number of predictors (single predictor, multivariate, correlated predictor errors); the foot points are no longer found by a separate inner optimizer, and orthogonality of the foot points is now the stationarity condition of the joint problem
* the Jacobian of the joint problem is now supplied to the solver exactly, using symbolic derivatives of the model formula (deriv()) and falling back to central finite differences if the formula cannot be differentiated (e.g. calls to user-defined functions). MINPACK's internal forward-difference Jacobian is unsuitable here, because its step is proportional to the size of each unknown and vanishes for small foot-point corrections, which stalled foot points at their observations (non-orthogonal points, e.g. for the last point of the odrpack guide example 2.C.i)
* control$outer_max is now the total budget of Levenberg-Marquardt iterations (default 5000). As minpack.lm silently caps a single nls.lm call at 1024 iterations, the solver is restarted from its last iterate until it converges, a restart brings no further improvement, or the budget is used up
* new control elements `ftol` and `ptol` (convergence tolerances, default 1e-10 instead of sqrt(.Machine$double.eps), which left points with very small residuals slightly non-orthogonal), `diag` (optional parameter scaling, MINPACK's adaptive scaling is now the default) and `lower`/`upper` as alternatives to the arguments of the same name
* `extend` and `window` now default to NULL, i.e. foot points are unbounded as in ODRPACK. They remain available as optional foot-point bounds for single-predictor models (`window` only together with `extend`)
* the covariance matrix of the parameters is now the ODRPACK (Gauss-Newton) covariance, i.e. the parameter block of the inverse of J'J of the joint problem, evaluated with effective-variance weights w_i = 1/(1/Qyy_i + g_i' * Qx_i^-1 * g_i) with g_i = df/dx at the foot point. This is valid for any number of predictors (the previous multivariate version ignored the predictor error in the standard errors)
* the NLS warm start now substitutes fixed parameters into the model formula instead of fitting them freely and resetting them afterwards, and no longer fails silently when `weights` are supplied (the weights were not found by nlsLM, so the raw starting values were used instead)
* onls() now warns if some free parameters have no measurable influence on the fitted model (e.g. a sigmoid that has collapsed to a constant because its exponent underflows), a degenerate stationary point at which the solver may report convergence
* the model formula is now evaluated in its own environment, so that user-defined functions and constants used in the model are found; the model is evaluated once at the starting values and errors are reported instead of being swallowed
* extended the core algorithm from single-predictor to n-dimensional (multivariate) models
* added support for weighted/heteroscedastic orthogonal regression via new `sigma_x`, `sigma_y` and `known_sigma` arguments (scalar, per-predictor vector, full p x p covariance matrix, or observation-specific standard deviations)
* check_o() now automatically switches from the tangent-angle criterion to a KKT-residual criterion for weighted fits, and checks multivariate models one predictor axis at a time
* check_o() no longer returns NaN if a foot point coincides with its observation (dx = 0, e.g. where the model slope is zero): the angle is now computed as atan2(|m*dx - dy|, |dx + m*dy|), which is algebraically identical to the old formula but avoids the division by dx. Observations lying on the fitted curve are reported with alpha = NA and Ortho = TRUE
* check_o() gained the arguments `tol_deg` (default 0.05 degrees) and `tol_kkt` (default 0.001) for the orthogonality tolerances, and now computes its slopes itself (central differences, evaluated in the environment of the model formula) instead of relying on internal helper functions
* plot.onls() gained an interactive 3D rgl-based plot for two-predictor models and a grid of partial-dependence panels (selectable via new `panel` argument) for more than two predictors
* plot.onls() gained an `asp` argument to toggle the 1:1 axis-ratio enforcement used to visually verify orthogonality, and now auto-derives a matching range when only xlim or ylim is supplied
* logLik_o() now includes a Jacobian/normalizing-constant correction for response/predictor precisions, so AIC()/BIC() from it are valid across onls() fits with different weighting schemes
* confint.onls() now explicitly conducts bootstrap confidence interval estimation on the orthogonal residuals. 
* fixed a diagonal p x p `sigma_x` matrix giving wrong (partly zero) predictor precisions
* fixed `weights`, `sigma_x` and `sigma_y` given as vectors not following `subset` and `na.action` (only a `sigma_x` matrix did)
* fixed residuals_o() returning values in sorted-predictor order instead of original observation order, which also misaligned NA reinsertion via na.action
* fixed summary.onls() reporting Std. Error/t value/Pr(>|t|) as 0/Inf/a wrong p-value for fixed parameters instead of NA, and computing cov.unscaled from the vertical- instead of the orthogonal-residual variance
* fixed check_o() using the wrong sign in its KKT-residual computation for multivariate weighted models, and now correctly realigns precisions to sorted-predictor order for single-predictor models
* fixed vcov.onls() rescaling the covariance matrix by the vertical- instead of the orthogonal-residual variance
* fixed logLik.onls() scaling the weight-correction term by n/2 instead of 0.5, affecting AIC()/BIC() for any fit with non-constant weights
* fixed predict.onls() erroring on se.fit/prediction intervals for models with fixed parameters
* fixed profile.onls() silently treating originally-fixed parameters as free while profiling other parameters
* deviance_o(), residuals_o() and print.onls() now report general precision-weighted orthogonal distances rather than assuming plain Euclidean distance
* removed the `...`, `jac` and `verbose` arguments from onls() (`trace` already covered progress reporting; the Jacobian is now always supplied internally)
* `data` is no longer optional in onls() (previously defaulted to parent.frame())
* onls() now stops with an informative error if the response is not a plain variable name, if all parameters are fixed, or if `lower`/`upper` do not have one value per parameter
* clarified that multivariate formulas must use explicit classical notation (y ~ b0 + b1*x1 + b2*x2) rather than R shorthand (y ~ x1 + x2)
* rewrote onls.Rd to describe the one-pass algorithm with its mathematics (objective, joint residual vector and Jacobian structure, stationarity conditions, covariance), updated check_o.Rd, and corrected the value of the bounded fit in example 6 (1.4376 / 0.9)
* added new reference-check examples to onls() (Total Least Squares closed-form comparison, Pearson/York weighted regression, Daeron & Vermeesch 2024, control settings) and new 3D/multivariate examples to plot.onls()

onls 0.1-4 (01-09-2025)

* fixed a URL error to retain CRAN compatibility

onls 0.1-3 (18-08-2025)

* fixed some small bugs to retain CRAN compatibility

onls 0.1-2 (31-10-2022)

* changed some URLs
* changed package maintainer
* improved grammar and style in manual and vignette
* removed <<- from onls.R
* changed NIST.R as to not print to console
* made the package more CRAN complient

onls 0.1-1 (07-09-2015)

* improved grammar and style in manual and vignette.
* transferred .Rmw vignette into 'vignette' folder for auto-generation    

onls 0.1-0 (26-07-2014)

* first package version
