########################
version 1.1 (2016-01-24)
########################

  Initial release, the sun shines and sommer has arrived!

########################
version 1.2 (2016-02-07)
########################

  + 'fdr' function bug was fixed
  + addition of the 'randef' function
  + addition of the converter 'atcg1234' function
  + names in the blup's or random effects added
  + zero-boundary constraint added to Average Information algorithm
     - it finds which var.comps are pushed to zero constantly
     - recalculates variance components removing such components
     - fix those values and calculates the most likely value for
       the problematic var.comp
  + now 'mmer2' can handle missing data in explanatory variables as lmer
  + now summary of 'mmer2' has names in the variance components
  + A.mat, D.mat and E.mat supported for polyploids
  + mmer can run GWAS for polyploid organisms
     - the models implemented are the same than Rosyara (2016):
     - "additive","1-dom-alt","1-dom-ref","2-dom-alt","2-dom-ref"
  + eigen decomposition to accelarate genomic prediction based on Lee (2015)
    has been added in the argument 'MTG2' of the AI, mmer and mmer2 algorithm

########################
version 1.3 (2016-03-09)
########################

  + The 'bag' function for bagging-GBLUP from Abdollahi-Arpanahi et al. (2015)
    has been added:
        - The function takes a model fitted and creates a bag matrix with
          the top markers (most significant) and creates a design matrix
          to be used as fixed effects in the GBLUP model to increase
          prediction accuracy.
  + 'bag' function has been equiped with stepwise selection to make sure that markers
    selected by "clustering" or "maximum" p.values methods provide at least a minimum
    increase in the prediction accuracy.
  + The Fisher Information matrix can be returned from the mmer function when
    the AI is used (default) but the argument 'Fishers' needs to be set to TRUE.
  + The bug for the AI algorithm when one var.comp and K and Z are diagonal has been
    fixed by changing to EMMA in this naive situation.
  + AI algorithm has been debuged to return the most likely variance components when the
    likelihood takes values around the maximum in a zig zag pattern. Just takes the value
    where the ML was found. When the likelihod follows a scale and dropping pattern
    the program will do the same. A warning message is emmitted.
  + GWAS modality of 'mmer' now adds the names of the markers of each score to keep track
    the value for each marker.

########################
version 1.4 (2016-03-24)
########################

  + The AI algorithm will take 5 EM steps if after 10 iterations (AI) the likelihood drops
    suddenly, indicating that initial values were too far from real values causing a bad
    behavior of the likelihood. The EM steps aim to provide initials values for those
    problematic variance components. ONLY mmer2!!!!!!
  + The AI likelihood behaving in a zig zag pattern is detected only after 10 iterations
    and we opted for returning the ML estimators.
  + Minor bugs have been fixed (names in random terms). In addition, ordering random effects
    based on their degrees of freedom has been implemented to provide more stability to the
    AI algorithm.

########################
version 1.5 (2016-05-03)
########################

  + Problem of "not-full rank X matrix" has been solved by reducing the X matrix column by
    column until is solved.
  + New alleles for deletions "-" has been added to the atcg1234 function
  + More examples have been added to the software, including maize Technow et al. (2015) data
    and rice Zhao et al. (2011)
  + We have opted for using the "EM" algorithm as default for the mmer2 function given the fact
    that most users using the mmer2 function would not use covariance structures. The mmer keeps
    using the AI algorithm as default method.
  + Implementation of the initial version of TP.prep function designed to select the best
    training population for genomic selection.

########################
version 1.6 (2016-05-19)
########################

  + Addition of the Newton-Raphson algorithm to sommer.
  + Vignettes with several examples available by typping vignette("sommer").
  + Bug in Principal Components within the function TP.prep has been fixed.
  + Example of GWAS for single cross hybrids has been added to the Technow data examples.
  + Bug on map.plot function for genetic maps fixed when having a column "Locus" as factor
  + Addition of winter bean data "FDdata" to exemplify the full diallel design calculation
  + mmer2 function now is fully functional as mmer but using data frames. All parameters
    can be manipulated as mmer except for the hdm function for half diallels.
  + The program will subset W if some markers are not polymorphic
  + beeper has beed added to the program to choose what sound to display when the program
    is don doing the calculations.
  + Now genetic maps can be added for a better display of the manhattan plots.
  + QQ plots are now displayed when doing GWAS for different models.

########################
version 1.7 (2016-05-23)
########################

  + Now the AI algorithm checks the situation when there are random effects with n < p
    (less observations than parameters to estimate) to take the right values in AI
    when extreme variance components values are found.
  + Bug in eigen decomposition of AI has been fixed by providing smaller starting values for
    the variance components
  + In addition, when n < p is present the initial values of the variance components are
    reduced significantly to provide better starting values.
  + New bug in the FDR function fixed. Now GWAS with or without map displais the FDR line.

########################
version 1.8 (2016-06-06)
########################

  + Minor changes, mainly in plotting defaults, warning messages, etc.
  + Addition of the phase.F1 function to create parental maps for F1(CP) crosses.
  + Addition of F1geno data accompaning phase.F1 function.
  + Argument MTG2 has been renamed as EIGEND.
  + EIGEND feature added to EMMA algorithm.
  + Bugs for errors when V(e) tends to zero have been fixed for all algorithms, warning
    messages have changed.
  + The function for finding biggest peaks "maxi.qtl" has been added.
  + Multiple responses can be fit using parallelization with the 'n.cores' argument

########################
version 1.9 (2016-07-01)
########################

  + lmerHELP has been slightly modified. Now only uses lmer initial values when required
    and only used for non-square random effects.
  + AI has been modified to dimminish updates values when they scalate to quickly.
  + zig zag detection in AI2 has been implemented, now even when boundary constraint is
    applied there should not be any problem.
  + Now 'bag' function respects if a map is available and does the bagging based on the map.
  + 'manhattan' function added to the package.
  + 'eigenGWAS' function included based on Chen et al. (2016) in Heredity.
  + Bug in AI failing when only one variance compnent was present and was zero, making the
    program collapse when trying to recalculate has been fixed.
  + Bugs in 'TP.prep' function fixed.
  + Now 'EMMA' algorithm handles multiple random effects.
  + The first multivariate algorithm "EMMAM" has been added to sommer (a single random effect).
  + Vignettes have examples for multi-response model and parallel univariate models.
  + 'bag' function has changed name to 'hits' to avoid confusion of selecting top 10 hits with
    Bootstrap aggregation (bagging).

########################
version 2.0 (2016-08-01)
########################

  + Multivariate algorithms "M-EMMA", "M-AI", "M-NR" have been implemented.
  + Minor bugs fixed.
  + Multivariate GWAS inplemented for any of the methods.
  + The function "mmer2" can run multivariate and parallel models as well.
  + Bug for parallel models fixed.
  + Eigen decomposition of additive relationship matrix implemented for
    multivariate models in "NR" and "AI".

########################
version 2.1 (2016-09-01)
########################

  + Standard errors and Z ratios added to the summary of mixed models.
  + Addition of functionality for residual structures, only AR1, CS, and ARMA
    supported. The most flexible functions "MNR" and "NRR" are not able to deal
    with missing data so you should be careful because the program imputes with
    the mean.
  + atcg1234 has a new functonality to get a presence/absence matrix for each
    allele at each marker.

########################
version 2.2 (2016-10-01)
########################

  + The LD.decay function has been added to sommer which takes a marker matrix
    and a map and calculates the LD decay.
  + In addition, the LD.decay function has the "unlinked" and "gamma" arguments
    that estimate the interchromosomal threshold for the gamma percentile to
    determine the real LD decay combined with the loess regression.
  + Bugs in GWAS models for univariate and multivariate forms when missing data
    exist have been fixed.

########################
version 2.3 (2016-11-01)
########################

  + The nearest neighbor function 'nna' has been added to adjust for neighbouring
    plots based on Lado et al. (2013).
  + The dataset 'ExpDesigns' with several datasets to teach users how to analyze
    certain experimental designs relevant to plant breeding has been added.
  + A fatal bug in mmer2 and mmer of imputing data when missing data existed has
    been corrected

########################
version 2.4 (2016-12-01)
########################

  + The IMP argument has been added to the function to allow the user to decide if the
    Y matrix should be imputed or get rid of missing values when estimating the variance
    components in the multivariate mixed models (only). The default is FALSE which means
    that the missing values are removed.
  + The gryphon dataset was included in the package to provide some help to the users
    that want tosee how pedigree data is used

########################
version 2.5 (2017-01-01)
########################

  + Bugs in fdr function fixed
  + Q+K model univariate fixed
  + Q+K model multivariate enabled
  + EIGEND feature added to NR method

########################
version 2.6 (2017-03-01)
########################

  + 'imp' argument added to the atcg1234() function to allow users to avoid imputation
  + Covariance structures for the residual component are not longer supported.

########################
version 2.7 (2017-05-01)
########################

  + at(.), diag(.), and(.), g(.) functions added to be used in mmer2.
  + NR algorithm updated to deal with multiple variance components equal to zero

########################
version 2.8 (2017-06-01)
########################

  + addition of the pin function to the package
  + issue with Year:g(id) type of arguments solved
  + more efficient EM algorithm; covariance matrices are inverted only once
  + EM fixed to don't calculate statistics at the end by direct inversion
  + rcov argument with options 'units' and 'at(.):units' enabled
  + AI MME-based added to sommer
  + new argument DI=TRUE/FALSE for deciding between MME-based and
    Direct inversion algorithms

########################
version 2.9 (2017-07-01)
########################

  + good versions of blocker and fill.design implemented
  + bug fixed in PEV for D-AI algorithm
  + minor bug fixes from summary functions in parallel models
  + going back to zero constraint where a model is fitted again witout the random
    effect close to the boundary.

########################
version 3.0 (2017-09-01)
########################

  + D.mat function now provides 3 different calculation methods
  + cleaner version of NR and AI algorithms
  + fixed EIGEND feature
  + AI algorithms now takes EM steps for better initial var.comp values
  + MME-base algorithms removed from mmer functions
  + remove of several function
  + move of GWAS functionality to a different function called GWAS and GWAS2
  + change of examples to show the flexibility
  + implementation of us(trait) and diag(trait) functionalities in mmer2

########################
version 3.1 (2017-11-01)
########################

  + pin function improved to don't work with scaled parameters
  + dominance relationship matrices corrected
  + general maintenance to several algorithms for minor bugs
  + constrained parameters in multivariate models not included in the summary as
    expected
  + better documentation about multivariate models

########################
version 3.2 (2018-01-01)
########################

  + variogram and plot.variogram functions enabled to visualize the residuals or the
    spatial model fitted (in case of 2D spline models)

########################
version 3.3 (2018-03-01)
########################

  + spl2D function to fit 2 dimensional spline model has been added and documented
  + spatPlots function added to visualize the fitted values and residuals
  + now the mmer and mmer2 function can use weights

########################
version 3.4 (2018-06-01)
########################

  + and() function has been replaced by the overlay() function effectively.
  + spl2D() can be used directly in the formula solver mmer2().

########################
version 3.5 (2018-07-01)
########################

  + h2.fun returns now corrected heritabilities for Oakey (2006) and Cullis (2006) formulas.

########################
version 3.6 (2018-09-01)
########################

  + some documentation changes
  + almost finished c++ version, just wait for it.

########################
version 3.7 (2018-09-01)
########################

  + C++ (Armadillo library) implementation done.
  + now the vs() function is the main function to create variance models
  + ds(),us(),cs(),at() for complex variance models added
  + constraint matrices can be easily indicated with unsm(), fixm(), uncm()
    in the Gtc argument of vs()
  + fatal bug in the multivariate BLUPs has been fixed
  + implementation of predict() function for lsmeans
  + multivariate GWAS implemented again, fixed and ready to rock
  + random regression implemented
  + anova() function enabled for single model to get sum of squares and MS
  + EIGEN functionality not available anymore, same with EMMA called from mmer
  + mmer2 and GWAS2 have been deprecated, now mmer and GWAS functions can do all

########################
version 3.8 (2019-01-01)
########################

  + stable version of c++ implementation, all bugs reported by users in 3.7 have been fixed
  + implement overlayed and leg in fixed effects

########################
version 3.9 (2019-04-01)
########################

  + documentation on the changes from old to new versions of sommer added to vignettes among better
    structure of the documentation
  + atcg1234() function now can take a marker matrix and a matrix with reference alleles and do the
    conversion with the customized reference alleles.
  + bug in the fixed effects depending of the order has been fixed. Is due to a strange behavior of the
    model.matrix() function.
  + Now Gtc and Gt arguments can take a list of matrices to apply different constraints when using the
    ds() and us() functions, although the cs() can provide the same results.
  + Now unsm(), uncm(), fixm(), fcm() have the rep argument to repeat the constraint matrix multiple times
  + Now DF for

########################
version 4.0 (2019-07-01)
########################

 + No updates for now other than documentation
 + Bug in the GWAS function for P3D=TRUE/FALSE fixed

########################
version 4.1 (2020-06-01)
########################

 + Bug for removing the intercept in models has been fixed
 + Better documentation for GxE models
 + Addition if the unsBLUP function to extract right BLUPs for unstructured models
 + Change of the vs() argument from Gt to Gti to specify initial values

########################
version 4.1.1 (2020-10-01)
########################

 + Redesign of the predict function to predict full grids
 + Better structure of the documentation
 + Bug in predict for models that include spl2D() or overlay() functions fixed

########################
version 4.1.2 (2021-02-01)
########################

 + Addition of the fitted() function to return Xb + Zu.1 + ... + Zu.n instead of just Xb
 + Bugs reported for the predict() function addressed
 + The function residuals() now returns e = y - Xb - Zu instead of e = y - Xb
 + A.mat, D.mat and E.mat functions are now C++ Armadillo implementations
 + Addition of the H.mat function in C++ Armadillo
 + When providing and Gu argument in the vs() function, missing levels in the data vector are added
   to provide all BLUPs
 + The predict function now provides the correct standard errors
 + GWAS C++ implementation
 + Progress bar for GWAS function added

########################
version 4.1.3 (2021-04-01)
########################

 + Addition of the gvs() function to fit indirect genetic effects and other competition models
 + Addition of an example of indirect genetic effects to vignettes.
 + Bug in predict() when missing data in covariates (x) exist was fixed.

########################
version 4.1.4 (2021-07-01)
########################

 + Addition of the stepweight argument to mmer() to control the magnitude of the update of the AI and NR
   when the information matrix goes out of boundaries too quickly.
 + Addition of emupdate argument to mmer() to control the use of expectation-maximization updates as
   alternative to the use of second derivative methods (NR, and AI).
 + Addition of the percChange output to the mmer() function to check the % change in the variance comp.
 + Argument buildGu for the vs() function added to avoid rrBLUP models to become slow. This implied modifying
   the MNR.cpp function to understand when the matrix multiplications with K (covariance matrix) should be
   avoided.

########################
version 4.1.5 (2021-11-01)
########################

 + Addition of the spl2Db and spl2Dmats() functions to fit all models capable to fit by the SpATS package.
 + Rename of spl2D function to spl2Da plus change of names for the arguments.
 + Fixes to a*b fixed model issues and predict improvements thanks to Sam Rogers and Julian Taylor.
 + GWAS by GBLUP examples added to the QG vignette. In addition new vignette for spatial modeling.

########################
version 4.1.7 (2022-07-01)
########################

 + Addition of the W argument to allows users to input a weights matrix instead of only a vector.
 + Henderson mixed model equation version of the AI algorithm based on Jensen, Madsen and Thompson (1997) finally availble with the mmec function. Please try it.

########################
version 4.2.0.1 (2023-01-01)
########################
 + P3D argument enabled again
 + Bug in mmec() for using weights with simple residual structures has been fixed. Now stage-wise analysis should work well. The problem still exist in the mmer() function though, where variance components are restricted in a different way.

########################
version 4.2.1 (2023-04-01)
########################
 + Predict function for mmec() now works well and returns the same results than asreml.
 + The predict function for mmer() on the other hand has been modified to calculate C = t(W) Vi W and Ci = solve(C), where W=[X Z] but when Vi is singular it makes difficult to come up with standard error identical to the mme-based methods.
 + The overlay() function now works in the fixed part of the formula.
 + The predict function now works with models using the overlay function.
 + A bug in the atcg1234 function that was causing failure when few markers were used has been fixed.
 + Weights can be used now in the GWAS function.
 + The emWeight argument now uses a decreasing logarithmic scale to assign weights to each iteration of the AI algorithm to guarantee convergence as best as we can.
 + Bug in the spl2Dc function in mmec has been fixed

########################
version 4.3.0 (2023-05-11)
########################
 + Bug in the rrc() covariance structure function when using a relationship matrix has been solved.
 + The new function redmm() function to create a reduced model matrix to fit huge models has been added

########################
version 4.3.1 (2023-06-11)
########################
 + There was a bug where the A matrix was not being reordered when provided in the mmec() function. It has been fixed.
 + A function to calculate effective population size based on allele coverage has been added.

########################
version 4.3.2 (2023-08-01)
########################
 + Bug in isc function fixed when matrices have a single column.
 + redmm() function improved calling the RSpectra package
 + now the mmec() function can use the different covariance structures in the fixed effect part so overlay and at models can be fitted as fixed.
 + Bug in atc() function fixed when a single level was used.
 + The covc() function has been added to calculate covariance between random effects as long as the incidence matrices of the random effects have the same levels.

########################
version 4.3.4 (2023-04-01)
########################
 + thetaC and theta arguments in special functions in mmec are better documented and have some examples
 + the atcg1234backTransform function has been added to bring back from numbers to letters.
 
########################
version 4.3.5 (2024-08-01)
########################
 + we have added the tps function to produce tensor product splines from the TPSbits package in a more straightforward way to produce incidence matrices that can be used in the mmec solver more easily.
 + keep improving the documentation of parameters and examples in different functions
 + rrc function for reduced rank models changed the specification and a good example is provided.
 + mmec function now ensures that data is transformed to as.data.frame() to avoid issues with tibble objects as reported in github.
 
########################
version 4.3.7 (2025-02-01)
########################
 + enhancing documentation

########################
version 4.4.1 (2025-04-04)
########################
 + unified mme and di algorithms unified in mmes()
 + removed the first two vignettes to make sommer lighter.
 
########################
version 4.4.2 (2025-06-04)
########################
 + now the direct inversion algorithm takes as input the sparse matrices to avoid overloading R memory with dense matrices. We transform them internally to dense matrices to make operations.
 + bug in getPEV fixed. Now the argument works again.

########################
version 4.4.3 (2025-07-29)
########################
 + updated version with better documentation for multi-trait models
 
########################
version 4.4.4 (2025-12-31)
########################
 + MNR function exported and documented for users wanting to specify their matrices directly

########################
version 4.4.81 (2026-09-30)
########################
 + ai_mme_sp2 function created to solve computation issues in the henderson algorithm. 
 + both function mmes and mmer are back and trying to split the direct inversion from henderson once again
 + now we can fit many different covariance structures using mmes 
 + more vignettes to explain covariance structures and theory behind
 + D.mat(nishio=TRUE) now implements Nishio and Satoh (2014), using the orthogonal
   dominance coding of Vitezica et al. (2013): -2q^2, 2pq and -2p^2, scaled by
   sum((2pq)^2). The earlier code centered heterozygosity by its observed column
   mean and divided by sum((2pq)^2), so D matrices (and E.mat with A#D or D#D)
   will differ from those of earlier versions. nishio=FALSE (Su et al. 2012) is unchanged.
 + D.mat and E.mat help pages corrected to describe the formulas that are used.
 + GLMM (PQL): any family with a valid link, including links with negative
   mu.eta (e.g. Gamma inverse). New binm(successes, failures/trials) builds a
   single proportion response whose trials are used as binomial prior weights.
   Pearson dispersion reported for free-dispersion families; negative binomial
   theta estimated with pqlControl=list(estimateTheta=TRUE). Working-model
   llik/AIC/BIC are NA and likelihood ratio tests are refused for PQL fits;
   residuals() gains type="pearson".
 + BUG FIX: theta_se of the Henderson engine was the inverse of the EM-blended
   update matrix on the -2logL scale (about half the true variance). It is now
   2*inverse(AI) over free parameters, matching the direct engine; standard
   errors in covParNative/covparams_mmes_se change accordingly.
 + New wald_mmes(): incremental and conditional (marginality-respecting) Wald
   tests with denominator df "none", "residual", "satterthwaite" or "kr"
   (Kenward-Roger). anova(object) with a single model returns it.
 + predict.mmes: new levels argument stored in a Dtable list column, sed/avsed,
   pairwise comparisons (all pairs or against a reference level) with p-value
   adjustment and residual/Satterthwaite/KR df for fixed-effect contrasts.
 + BUG FIX predict.mmes: "include and average" divided each cell by the total
   number of cells of the term instead of the cells present for the prediction
   row (e.g. 1/12 instead of 1/3 for V:N); results now match emmeans.
   Prediction labels no longer lose characters matching the classify name.
 + mmes objects store call and inputArgs so model inputs can be rebuilt.
 + New metricm(): exponential, power, gaussian, spherical and circular
   distance-based correlations, isotropic or anisotropic (product per axis or
   geometric angle/ratio), with analytic derivatives (ASReml exp/iexp/aexp/
   gau/igau/agau/isp/cir). maternm() gains geometric anisotropy.
 + New strm(): covariance structure spanning several random terms (ASReml
   str), sigma2 * K_terms x K_inner x A, with any constructor over the terms
   and levels aligned on the union or on Gu. covm() is its two-term case.
 + New vcc argument of mmes(): equality and scaling constraints between
   covariance parameters (ASReml vcc), identified through the new vcParams
   table (returnParam=TRUE and fitted objects). Implemented in both engines
   by projecting each AI update on w = o + T phi; theta_se follows the
   constraints.
 + New meffects_mmes(): marker effects recovered from a GBLUP term (ASReml
   mef) reproducing A.mat() or a custom transform, with optional standard
   errors and back-solved GWAS tests.
 + Wald, residual df and PQL dispersion use sparse rank computations.
 + New stackTraits() converts wide trait columns to long format with an
   explicit record key. mmes() accepts that key in the residual structure to
   pair trait residuals without relying on row order.
 + New familym() supports mixed-family multivariate PQL models, including
   binomial/Poisson/negative-binomial traits with fixed residual dispersion and
   Gaussian traits with estimated residual variances. corgm() gains a
   heterogeneous-variance mode for residual trait covariance structures.
 + New covmatrix_mmes() reports term covariance and correlation matrices with
   delta-method standard errors. Default usm()/dsm() starting shapes use
   per-coordinate residual mean squares.
 + BUG FIX Henderson compact residual score used an incorrect zero
   log-determinant derivative, which could prevent complete multi-trait usm
   fits from converging. Weighted Kronecker residual derivatives now use the
   same inverse-Cholesky weighting as the effective residual covariance.
 + BUG FIX weighted Henderson fits with non-diagonal Kronecker residuals now
   match the direct engine for covariance parameter estimates.
 + New mmes(solveOnly=TRUE, covPar=) solves the mixed-model equations at known
   covariance parameters (genetic evaluation) with a matrix-free,
   block-Jacobi preconditioned conjugate gradient that never assembles C.
   covPar takes a fitted mmes object, a list of natural parameters or
   G0/R0 matrices, or a parameter table. 1M animals (pedigree A-inverse):
   ~10 s single trait, ~40 s three traits.
 + vsm() no longer builds a dense records x levels mask per covariance
   coordinate, and mmes() builds the fixed design from distinct covariate
   rows and pairs residual records without ave(); input preparation for
   large data is several times faster and uses far less memory.

########################
version 4.4.90 (2026-10-13)
########################

  Engine speed work for ai_mme_sp2 (details and measurements in
  src/ai-mme-sp2-acceleration-plan.md).
 + New mmes_control(experimental=list(profile="force")) records per-phase
   seconds and workload counters in fit$engineDiagnostics$profile. New
   profile_mme_workload() tabulates them; run_mme_benchmark() carries phase_*
   columns; benchmarks/profile_phase0.R profiles AR(p), SAR/CAR and dense-Gu fits.
 + New mmes_control(memoryPolicy="speed"|"limit"), default "speed": when the
   CHOLMOD selected-inverse caches do not fit denseMemoryMB the fit stops and
   reports the memory needed, instead of a silently slower fallback
   ("limit" keeps the old fallbacks).
 + The block-banded engine is no longer chosen automatically (up to 80x slower
   than generic CHOLMOD on AR1 fits); it remains available with
   experimental=list(blockBanded="force").
 + General dense block-Schur admission now compares predicted sparse-factor
   work (AMD symbolic fill) with dense-plan cost instead of block density.
   Coupled plans require a work/cost ratio of at least 0.05; uncoupled groups
   retain their solve/trace fast path. Invalid estimates do not admit Schur.
   An isolated low-fill random-graph fixture improved from 2.968 s to 0.641 s,
   with maximum relative output differences below 1.6e-9 across nine fixtures.
   Chain/latent/banded policies are unchanged by this guard. The CHOLMOD dense
   inverse-block cache still skips sparse groups. The provisional 64/1000
   fragmentation rule was removed after true on/off validation contradicted
   the earlier memory-budget proxy comparison.
 + The positive-definiteness check of single-factor covariance structures with
   >= 128 effects (SAR, CAR, large AR/Toeplitz) uses a shifted Cholesky test
   instead of an eigendecomposition; identical fits, 2-4x less time in the
   constraint step.
 + mmes() starting values for variance components no longer use a dense QR
   of the fixed-effects design (a 12000 x 1000 design spent about 90% of the
   wall time there with multithreaded BLAS); a sparse normal-equations solve
   is used for well-conditioned full-rank designs. Same converged estimates,
   about 6x faster end to end on that model.
 + New mmes_control(experimental=list(blockSchur="auto"|"off"|"force")) to
   disable (off) or bypass the density/group-size rules of (force) the block
   elimination engines; engineDiagnostics$profile gains traceFallbackSolves.
 + Opt-in experimental=list(traceBatch="force") batches exact LDLT random-score
   fallback solves in independent OpenMP panels, without the shared solve lock.
   A 1000-level AR2 fixture improves from 7.959 s to 6.084 s with four threads;
   likelihood, scores, AI, estimates and uncertainty match the old exact path.
   Peak RSS increases about 14.2 MiB. One-thread fits do not benefit, so the
   default stays off. Panel scratch obeys speed/limit memory policies.
   New profile counters: traceBatchSolves and traceBatchWorkspaceBytes.
 + Opt-in experimental=list(sparseFactors="force") uses exact sparse AR(p)/ANTE
   precision operands in derivative products and tensor applications, compact
   ANTE innovation rows, and sparse SAR/CAR covariance/derivative solves.
   On paired 1000-level fixtures, whole-fit medians improve from 6.287 s to
   4.770 s for SAR and 5.699 s to 4.168 s for CAR; peak RSS falls about
   114.9 MiB and 45.0 MiB, respectively. Likelihood, scores, AI, estimates
   and uncertainty agree within 5.23e-9 relative error across ten fixtures.
   AR/ANTE whole-fit gains are modest or unestablished, so defaults stay off.
   Dense compatibility and covariance-derivative storage remain; no derivative
   thresholding, PCG estimator change, or optimizer change is introduced.
 + LDLT selected inverse indexes dense columns by position instead of hashing
   (solver="ldlt" with a dense Gu, 300 individuals: 72 s to 16.5 s).

 + ai_mme_sp2 now delegates covariance evaluation, precision, derivatives and
   tensor operations to a shared covariance engine. Structured kernels use
   representation capabilities rather than model labels, including custom
  modified-Cholesky innovation callbacks. Saved built-in factors acquire
  missing capabilities at the R compilation boundary. Path and workspace policy remains
   in mme_paths.h; existing acceleration defaults and output names are retained.
 + BUG FIX complete-block ANTE residuals with sparseFactors="force" no longer
   reject valid covariance derivative indices. Multi-factor SAR/CAR scores now
   include the Kronecker log-determinant multiplicity. Sparse spatial evaluation
   caches obey their aggregate entry allowance without unsigned underflow.

 + rrm() now fits genuine reduced-rank covariance with all specific variances
   fixed at zero, matching the covariance family of ASReml rr(). The first
   loading is fixed to one to identify the overall scale supplied by vsm().
   Native latent-factor Henderson and direct-inversion fitting support the
   singular random covariance without an artificial diagonal remainder.
 + fam() gains variance="homogeneous" to estimate a common specific variance
   (the former rrm() covariance family); the default remains "heterogeneous".
   loadings_mmes() returns exactly zero specifics for RR and scales specific
   variances by sigma2, not sqrt(sigma2). RR regression scores recover latent
   BLUP coordinates; Bartlett scoring rejects zero-specific RR models.
 + Native latent RR fitting now uses loading-aware scores and average
   information instead of BFGS and numerical observed information. Covariance
   derivatives include changes in the latent design; the original
   covariance-family descriptors and sparse Henderson solves are reused.
   Parameter uncertainty uses final unblended AI, while stepWeight and
   emWeight control updates and the iteration monitor records accepted fits.
 + Native RR AI traces now share exact design and residual contractions across
   covariance parameters, avoiding repeated observation-space derivative and
   relationship solves. Conditional likelihood trials skip unused ordinary
   trace preparation. Verbose RR fits report their MME path, coefficient solver
   (including eligible block-Schur dispatch), and accepted iteration progress.
 + RR progress prints the coefficient engine before the table header and
   includes wall/CPU times, EMweight and actual accepted trust/line-search
   scales. Numeric values are retained in engineDiagnostics$iterationProgress;
   inapplicable boundary-restraint and notPD columns remain omitted.

 + BUG FIX native RR AI fitting now supports dsumm() spatial residuals.
   Covariances, analytic derivatives, and shared trace contractions operate
   on observed residual blocks rather than expanding the section-by-grid
   covariance. Existing section slicing and factor derivatives are reused.

 + Native RR AI now reuses the live LDLT/CHOLMOD or block-Schur coefficient
   solver and the likelihood's residual precision operator. Spatial fits use
   the shared grid engine and its factor-wise derivatives rather than a
   second dense observed-block precision. Random traces cancel relationship
   covariance actions through the latent Henderson identity. Verbose output
   identifies the shared residual engine; AI phase timings are retained.

 + Opt-in mmes_control(acceleration="hybrid-ai") reduces repeated Henderson
   sensitivity solves using Powell-damped BFGS between exact AI refreshes.
   Scores, likelihood acceptance, bounds and PD safeguards remain unchanged;
   intermediate stabilization uses the BFGS diagonal. Final InfMat and
   covariance-parameter uncertainty always use exact AI, including at the
   iteration limit. Convergence requires exact step and projected-score checks.
   Gaussian LDLT/CHOLMOD fits are supported; vcc ties explicitly retain exact
   AI. Other optimizers and the default behavior are unchanged. Diagnostics
   expose refreshes, resets, certifications and sensitivity RHS counts.

 + Experimental mmes_control(experimental=list(saddlePoint="force")) adds an
   exact augmented Henderson MME solver using pivoted Eigen sparse LU.
   Defaults and "auto" remain unchanged. Backend-independent exact panel
   contractions and normal-equation output assembly are reusable in
   mme_exact.h; the augmented factorization lives in mme_saddle.h.
   REML/ML retain exact scores, AI, covariance SE and coefficient uncertainty.
   Fitting avoids normal-equation cross-products; final C is assembled once
   for prediction and saved-object compatibility. Unsupported routes fail
   explicitly. Initial benchmarks do not establish a speed or memory gain;
   the solver remains force-only, with no automatic admission rule.

 + mmes() now blends genuine EM/GEM and AI proposal directions rather than
   information matrices. blendWeight replaces emWeight (a deprecated named
   alias), and blend="em" selects this method.
   Scalar, fully free diagonal and unstructured terms use closed-form
   posterior-moment updates when available; other structures and constraints
   use numerical increasing-Q generalized EM. Complete unstructured residual
   blocks also have closed-form M-steps. Mixtures use likelihood
   backtracking with an EM fallback. Full EM weights skip iterative AI,
   zero weights skip EM, and final inference always uses exact AI.
   Positive weights require deterministic Gaussian Henderson REML or ML without
   other accelerators, factor-score augmentation, PQL or saddle-point fitting.
   Unsupported modes default to pure AI when blendWeight is omitted;
   explicitly positive weights remain errors.
   Blended fits default to stepWeight=1; pure AI retains its previous damping.
   engineDiagnostics$blend reports work counts, Q gains and numerical
   M-step convergence. Exact moments obey the denseMemoryMB budget.
   Henderson weight factorizations now consistently preserve observation
   order across EM and AI. This fixes a previous correlated-weight ordering
   discrepancy; affected weighted estimates can differ from earlier versions.
   A scalar whole-number blendWeight greater than 1 now sets the exponential
   decay length from 1 to 0.03, capped at nIters, followed by constant 0.03.
   Fractional counts and vectors containing values above 1 are rejected.

 + blend="pxem" now selects scalar parameter-expanded generalized EM (PX-GEM)
   blended with AI. Each random term has a positive expansion scalar; joint
   numerical optimization preserves all covariance families, fixed values,
   scale floors and vcc constraints in original reduced coordinates.
   Existing EM/GEM updates seed the expanded M-step. Posterior moments include
   fixed/random and cross-random uncertainty, and reuse existing MME paths.
   PX proposals use observed-likelihood backtracking with ordinary EM/GEM
   fallback. Full weight skips iterative AI; zero weight skips EM/PX work.
   Final inference remains exact AI. Diagnostics expose expansion factors,
   additional Q gains, numerical convergence/evaluation counts and fallbacks.
   Additional moment/optimizer workspace obeys denseMemoryMB. This is scalar
   numerical PX-GEM, not full matrix expansion; fewer iterations do not
   guarantee lower runtime. The default remains blend="em".
   Reduced-rank EM/PX-GEM now cache exact loading-basis posterior moments
   once per E-step instead of repeating coefficient solves at every numerical
   M-step objective probe. This fixes severe first-iteration latency in
   models with many loading and residual parameters, without changing Q.
   EM/PX random moments now reuse exact Schur/inverse-block and selected-inverse
   contractions, requesting only diagonal moments for diagonal covariance
   updates. Scalar residual moments reuse a fixed-effect-aware trace identity.
   Scalar non-reduced-rank PX expansion moments reuse sparse coefficient
   contractions when memory permits. Generic bounded exact panel solves remain
   the fallback. New work counters expose contraction and fallback usage.
   AI directions inside EM/PX schedules now use native exact AI kernels for
   non-reduced-rank models instead of generic latent information calculations.
   Zero-weight iterations skip moment callbacks. Mixed iterations share native
   AI factorizations and inverse caches with posterior-moment work. Accepted
   exact-information trials supply subsequent zero-weight iterations and final
   inference without duplicate factorization. Reduced-rank AI retains its shared
   latent kernels; blend guards and constraints remain unchanged.

 + Omitted blendWeight now selects AI-first optimization with on-demand
   EM/GEM recovery (or scalar PX-GEM with blend="pxem") for eligible Gaussian
   Henderson REML or ML fits. Healthy AI iterations skip EM/PX moment and objective
   work. Repeated poor progress, failed AI line search, or unreliable
   information solves trigger one recovery step when the projected score
   remains non-small. Recovery must improve likelihood before AI resumes.
   Automatic convergence requires a scaled projected-score certificate.
   Failed recovery warns and returns the last accepted point without
   convergence. Explicit numeric blend schedules and unsupported-engine
   defaults retain their behavior. Recovery diagnostics record triggers and
   stop reasons; automatic refits retain blendWeight=NULL.

 + ML now supports explicit EM/PX-GEM schedules and default AI-first recovery.
   ML posterior moments use the random-only D inverse, with deterministic
   current GLS fixed effects during covariance and loading Q/PX probes.
   Likelihood trials reprofile fixed effects and accept against profile ML.
   Random, residual, weighted and reduced-rank moment paths share a separate
   posterior action; joint solves and final coefficient uncertainty are unchanged.
   CHOLMOD retains independent C/D symbolic caches across conditional fits.
   Diagnostics identify ML versus REML. Dense conditional-moment and expanded-Q
   tests cover both likelihood modes.

 + mmes() iteration output now shares one compact formatter across ordinary,
   hybrid and reduced-rank AI, scheduled EM/PX blends, AI-first recovery and
   direct inversion. Columns show accepted likelihood/gain, update method,
   recovery policy/trigger, trust and line-search scales, notPD repairs, and
   elapsed seconds. Redundant wall-clock/CPU columns are removed; detailed
   stored diagnostics remain unchanged. Recovery rows identify the actual
   method tried, including late line-search recovery and EM/PX fallbacks.

 + Unreliable automatic AI solves now try a scaled, ridge-regularized,
   trust-limited AI direction before computing EM/PXEM moments. Both retries
   require feasible observed-likelihood improvement; ordinary AI resumes
   after each accepted step. Explicit numeric blend schedules are unchanged.
   Progress reports AI-reg and diagnostics record regularized attempts,
   accepted steps, source iterations, and ridges.

 + Standalone sparse, saddle, and PCG native tests now check for their private
   source headers before compiling. They skip in installed-package check
   layouts without src/, while source-tree runs retain the native checks.

 + Direct inversion (henderson=FALSE) now supports the same blend/blendWeight
   controls as Henderson: true observation-space EM/GEM, scalar PX-GEM,
   mixed directions, scalar/vector schedules, and exponential tapers.
   Omitted weights select lazy AI-first recovery; explicit zero preserves
   AI-only fitting. All direct covariance families support ML and REML,
   including custom factors, constrained and ragged residual structures,
   and zero-specific reduced-rank latent terms. Numerical M-steps increase Q;
   they do not claim global maximization. ML profiles fixed effects in Q,
   and REML includes joint fixed-effect posterior uncertainty.
   Full EM skips AI update solves, while final inference remains clean AI.
   Direct EM workspace obeys denseMemoryMB. Diagnostics record Q work,
   actual methods, numerical M-step status, recovery, and likelihood reuse.
   Henderson algorithms are unchanged; the existing numerical Q minimizer
   is shared without changing its implementation.

 + Opt-in experimental=list(reducedAI="force") now packs independent vcc
   directions before Henderson sensitivity solves and assembles exact reduced
   AI without full sensitivity or AI buffers. Default AI-first recovery and
   positive EM/PXEM direction blending support Gaussian deterministic REML/ML,
   preserving fixed/tied coordinates, boundary restrictions, trust units,
   regularized AI recovery and exact score certification. Full covariance
   reporting retains a full exact final evaluation. covarianceInformation="none"
   avoids report-only full information while preserving coefficient uncertainty.
   Default remains off. Auto admission reports unsupported-path fallbacks;
   force rejects them. Legacy explicit blendWeight=0, acceleration, latent
   reduced rank, PCG, PQL, solve-only, saddle-point and factor-score augmentation
   are not supported. Full score traces and active local curvature-pair work
   remain unchanged; small full-reporting fits can be slower. Diagnostics expose
   actual/skipped sensitivity columns and reduced workspace sizes.

 + Opt-in experimental=list(commonScale="force") profiles the overall covariance
   scale analytically in eligible Gaussian deterministic Henderson REML/ML fits.
   Independent shape sensitivities incorporate vcc ties directly; exact profile
   scores and AI scale Schur elimination need no extra factorization per candidate.
   The dedicated optimizer uses trust caps, regularized AI, likelihood backtracking
   and exact final score certification, not EM recovery. Product-level variance
   scales must be free, fixed designs full column rank, and ties common-scale
   compatible. Fixed shapes and positive variance-ratio ties are supported.
   Only interior product variances above max(1e-8,tolParInv) are supported:
   initial floor violations error, trial violations backtrack and stalls warn
   without declaring convergence. Refit with commonScale="off" for boundary
   fitting or ordinary AI/EM recovery. Positive blendWeight, PCG, PQL, solve-only,
   acceleration, latent reduced rank, saddle-point and factor-score augmentation
   are excluded; omitted or zero blendWeight is accepted. Default remains off.
   Auto admission exposes static fallback reasons rather than hiding runtime
   failures. A final ordinary evaluation preserves original-coordinate outputs,
   coefficient uncertainty and requested full covariance scale/shape information.
   Information omission still performs that final evaluation. Wald/KR perturbed
   fixed-parameter evaluations bypass scale profiling. Common-scale diagnostics
   expose total work and certification; no general speedup is claimed.

 + Breaking numerical-control API: mmes() now accepts one modular control
   object through control=mmes_control(). Numerical arguments and the former
   engine/experimental/pcg containers are replaced by task, formulation,
   representation, solve, logdet, score, curvature, update, proposal,
   convergence, output, execution and diagnostics settings. Unknown names
   are errors; legacy public aliases are not accepted. Model formulas,
   families, covariance constraints and known covariance inputs remain
   model arguments. See sommer.mmes.options for migration examples.
   Default algorithms, schedules, tolerances, work budgets and output requests
   retain the previous numerical bundle. Aitken extrapolation belongs to
   proposal; hybrid AI belongs to curvature. Their combination remains
   unsupported. Shared typed native policy separates those responsibilities
   without changing numerical kernels. Fitted objects retain control and an
   executionPlan; detailed native diagnostics remain the eligibility evidence.

 + mmes_control() and mmes() check model-independent option combinations
   against one compatibility table. Supplied settings that conflict, or that
   the selected path would ignore (for example PCG estimator settings with a
   deterministic solver, Henderson settings with formulation="direct", or
   estimation modules with task$mode="solve"), stop with one message naming
   every offending setting and a fix. Restated defaults are accepted.
   Changing task$pql on a Gaussian identity fit is an error. Fit-time
   eligibility messages now use module option names (e.g. solve$method,
   output$inverse, representation$equations) instead of removed arguments.

 + Consolidated legacy mmer, its methods, covariance constructors and GWAS
   functions in R/DEPRECATED_mmer.R without changing their definitions.
   Removed duplicate current-path helper definitions.
 + Current mmes post-fit prediction, posterior variance and marker inference
   reuse operation-local coefficient factorizations. Marker SE calculations
   form relationship-marker products by block, reuse retained precision and
   preserve sparse supplied inverses.
 + Native covariance reporting now propagates cross-term uncertainty jointly;
   vpredict uses the full native covariance matrix. Unavailable uncertainty
   does not contaminate independent native quantities.
 + Reduced allocations in covariance-matrix SEs, leverage diagonals, spatial
   splines and APY assembly; reused covariance-factor compilation, shared
   control defaults and PQL warm-start compatibility checks.
 + Benchmark results now record the actual convergence flag, retain uniform
   rows on warmup/timed failures, compare final likelihoods and mark unchecked
   or unavailable-reference parity as NA.

 + Simplified the regression suite by removing overlapping family/solver
   smoke sweeps, detailed internal-counter and output-format checks, and unused
   reference helpers. Retained independent numerical references, public API
   validation and reported-bug regressions. Replaced spatial hash snapshots
   with direct tensor/design checks and reduced repeated threading stress fits.
   Future tests should cover a distinct public contract or numerical failure,
   rather than repeat the full cross-product of already tested options.

 + Package consolidation: separated covariance descriptor/family code and
   formula helpers without changing their contracts. Lifecycle hooks now
   live in R/zzz.R. Native utilities, relationships, post-fit inference,
   direct fitting and known-covariance solving have separate translation
   units; shared direct covariance evaluation and MME design assembly have
   explicit internal boundaries.
 + Consolidated public-to-private PCG option mappings and added ownership
   and mutable-control regression tests.
 + Added cross-platform package-check CI and a required source-tree native
   test job. Existing standalone native control checks now run automatically,
   with temporary caches and restoration of both global test options.
 + Reduced benchmarks to one seeded performance/parity runner and its usage
   guide; removed historical outputs, installed libraries and exploratory
   scripts. Old recorded measurements remain historical, not current files.
 + Added an executable README quick start and contributor/test guidance.
   MASS and crayon are implementation imports rather than attached
   dependencies; use qualified calls (e.g. MASS::mvrnorm) in user scripts.
   Matrix and enhancer remain attached for the established model interface.

 + Fixed Henderson random-effect vector u to use the original response scale,
   consistent with bu and uList. Added independent mixed-model equation checks
   across LDLT, CHOLMOD, PCG and direct fitting, and strengthened prediction
   fixtures to require converged fits away from zero-variance boundaries.

 + Gaussian REML fits now omit independent random covariance components whose
   designs lie entirely in the fixed-effect space. The conservative check
   handles free single-axis dsm/atm variances and whole free-scale terms with
   fixed shapes. Omitted components retain zero covariance and random effects
   in public outputs, with NA uncertainty, an explicit status and an audit
   table; numerical regularization is no longer mistaken for estimability.

 + Added opt-in Gaussian REML conditional FA updates through
   mmes_control(update=list(method="conditional-fa")). Sparse posterior
   moments are a separate mathematical block with per-fit symbolic reuse,
   refreshed numeric factors, cached sparse quadratic contractions and
   bounded inverse panels. Incomplete/repeated records and diagonal or FA
   residuals are supported; unsupported combinations fail explicitly.
 + Added fam(shared="group") for common genetic/residual loading axes and
   source-specific diagonal factor variances on the conditional path.
   Factorized outputs omit dense trait covariances and uncertainty matrices.
   This exact EM foundation is experimental, not a thousand-trait speed
   guarantee; convergence may be slow near variance boundaries.
 + Regression FA scores now use a rank-sized Woodbury solve instead of a
   dense trait covariance solve, including factorized conditional-fit output.

 + Separated latent FA equations from optimization: use
   representation=list(equations="conditional") with update$method="ai"
   (default) or "em". The earlier method="conditional-fa" remains EM.
   Added exact analytic REML scores and average-information contractions,
   joint shared-loading derivatives, reusable posterior-factor RHS solves,
   rank-revealing safeguarded AI steps and explicit EM recovery diagnostics.
   Physical-specific/allocation bounds use execution$inverseTol and report
   their active KKT coordinates, preventing near-deterministic latent systems
   from accepting numerically unreliable likelihood improvements.
   Conditional AI admits dense derivative/curvature workspace through the
   existing memory policy; ordinary AI and conditional EM defaults are unchanged.

######################
## PENDINGS AND PRIORITIES
+ Fix plot.mmes
+ check how dense schur engine is used when more than one random term with dense Gu is used