`AnimalBreedingTools.approximate_trace_of_inverse`

— Methodtr = approximate*trace*of_inverse(A)

`AnimalBreedingTools.bending`

— Method```
bentB = bending(B [, W], gamma=nothing, tol=eps, verbose=false)
bentB = bending(B, gamma=nothing, tol=eps, verbose=false)
```

Apply a traditional "bending" method to a symmetric matrix `B`

(between-group-variance) using another symmetric matrix `W`

(within-group-variance), as described by Hayes and Hill (1981). The resulting matrix is `bentB = (1-gamma)*B + gamma*v*W`

, where `v`

is the average eigenvalues of `B`

and `gamma`

is a constant between 0 and 1, which can be provided by the user or generated by this function. When `W`

is omitted, the identity matrix will be used instead.

If the user supplies `gamma`

, the function will simply apply the bending formula and return `bentB`

without checking whether `bentB`

is positive defnite or not. Otherwise (`gamma=nothing`

), this function generates `gamma=1/1000`

, applies bending, and returns `bentB`

if it is positive definite as `minimum(eigen(bentB).values)>tol`

(the default `tol`

is the machine epsilon). If not, the function tries `gamma=2/1000`

, `gamma=3/1000`

, ..., and `gamma=1000/1000`

until `bentB`

becomes positive definite.

With `verbose=true`

, this function prints `v`

(average eigenvalue of `B`

) and `gamma`

that has been used for bending (default: `verbose=false`

).

`AnimalBreedingTools.bending2`

— Method`bentV = bending2(V [, W], corr=false, mineig=eps, tol=eps, maxiter=10000, force=false, verbose=false)`

Apply an alternative "bending" method to a variance-covariance matrix `V`

using a weighting matrix `W`

, as described by Jorjani et al. (2003). When omitting `W`

, `J`

(matrix of ones) will be used, implying the same weights for all elements. Jorjani et al. (2003) suggested `W`

as the reciprocal of number of individuals used in the estimation of covariances.

The function calculates the eigenvalues of `V`

, and an eigenvalue is replaced with `mineig`

if it is smaller than `mineig`

. The covariance matrix is reconstructed with modified eigenvalues, then tests whether it is positive definite by checking the minimum eigenvalue larger than `tol`

. The function repeats the process until the resulting matrix becomes positive definite, and the number of maximum iterations is defind by `maxiter`

.

See some options below.

- With
`corr=true`

, this function assumes the input is a correlation matrix; otherwise, the input is expected to be a variance-covariance matrix. - With
`force=true`

, this function iterates the formula`maxiter`

times regardless of whether it has converged or not. - With
`verbose=true`

, this function shows the number of iterations required when`force=false`

.

`AnimalBreedingTools.chol_symtrid`

— Method` dv,ev = chol_symtrid(T::SymTridiagonal)`

Calculate the Cholesky factor of a symmetric tridiagonal matrix. Note that the factor is a lower tridiagonal matrix, and the output is a couple of vectors. The factor is stored in `dv`

for the diagonal matrix (D) and in `ev`

for the off-diagonal elements (L).

`AnimalBreedingTools.directsum`

— Method` M = directsum(X...)`

Returns the direct sum of supplied matrices.

`AnimalBreedingTools.disp`

— Method```
disp(x)
disp(x...)
disp(io::IO,x)
disp(io::IO,x...)
disp(m::Array; format)
disp(io::IO, m::Array; format)
```

Write an object `x`

(or objects `x...`

) to `io`

(or `stdout`

if `io`

is not given). The output is the same as `println`

and `Base.print_matrix`

, but the element type will not be shown. If an array is given, you can specify the C-style format with `format`

.

```
A = [1 2; 3 4]
b = [9,10]
disp(A)
disp(b)
disp(A,b)
```

`AnimalBreedingTools.get_design_matrix`

— Method```
Z = get_design_matrix(x, nested=[]; maxlev=0, cov=false)
Z = get_design_matrix(x, y)
Z = get_design_matrix(X; maxlev=[], cov=[])
```

Returns a design matrix accoding to class variables or covariates in a vector `x`

or matrix `X`

. For a vector `x`

, optional argument `dested`

defines `x`

as the nested covariate within a class. Option `maxval`

specifies the number of levels in this effect. Omitting `maxval`

, the function uses the maximum integer as the maximum level in `x`

or `X`

. Option `cov`

specifies the input is covariate and passes it through the output. The matrix `X`

doesn't support nested covariates.

If two vectors `x`

and `y`

are provided, the function returns the design matrix for the intercation between the two effects; the main effects are not included in the result. Also, in this case, the function does not accept any additional options.

```
julia> x = [1,2,1,1,2]
julia> X = get_design_matrix(x)
5×2 Array{Float64,2}:
1.0 0.0
0.0 1.0
1.0 0.0
1.0 0.0
0.0 1.0
```

`AnimalBreedingTools.get_interaction_matrix`

— Method`X = get_interaction_matrix(X1, X2)`

Returns a design matrix for the interaction between effects. The design matrices for the two main effects (`X1`

and `X2`

) should be provided. Note that this function will not check if the supplied matrices are really the design matrices - it just creats all the combination of column products among the input matrices.

`AnimalBreedingTools.hrv_class`

— Method`vare = hrv_class(ve,rangeset)`

Generate a vector containing residual variances from a vector of heterogeneous residual variances `ve`

and a set of ranges `rangeset`

.

```
julia> vare = hrv_class([10.0,20.0],[3:5,6:9])
7-element Vector{Float64}:
10.0
10.0
10.0
20.0
20.0
20.0
20.0
```

`AnimalBreedingTools.ldlt_symtrid`

— Method`dv,ev = ldlt_symtrid(T::SymTridiagonal)`

Calculate the LDLT factor of a symmetric tridiagonal matrix. Note that the factor is a lower tridiagonal matrix, and the output is a couple of vectors. The factor is stored in `dv`

for the diagonal matrix (D) and in `ev`

for the off-diagonal elements (L).

`AnimalBreedingTools.legendre_polynomials`

— Method`L = legendre_polynomials(n; normalize=true)`

Calculates an upper triangular matrix with the `n`

-th order Legendre coefficients. By default, the coefficients will be normalized as `normalize=true`

.

`AnimalBreedingTools.load_data_set`

— Function`data = load_data_set(name)`

Load a small data set for numerical testing.

`AnimalBreedingTools.normalized_legendre`

— Method`v = normalized_legendre(x,n)`

Calculates a single value of the `n`

-th order Legendre polynomial of `x`

. It calls `Jacobi.legendre`

.

`AnimalBreedingTools.normalized_legendre_matrix`

— Method`M = normalized_legendre_matrix(t,n; tmin=minimum(t),tmax=maximum(t))`

Calculates a matrix with the `n`

-th order Legendre coefficients on time points `t`

, a vector or range. The earliest and latest (minimum and maximum) time points will be taken from the edges of `t`

. You can provide these time points, `tmin`

and `tmax`

, as options.

```
# up to 3rd order coefficients for 5 to 12
julia> M = normalized_legendre_matrix(5:12,3)
# same expression
julia> M = normalized_legendre_matrix([5,6,7,8,9,10,11,12],3)
# as above with time points ranging 1 and 15
julia> M = normalized_legendre_matrix(5:12,3, tmin=1, tmax=15)
```

`AnimalBreedingTools.r2v`

— Method`cov = r2v(corr, stdev)`

Convert a correlation matrix `cov`

and the vector of standard deviation `stdev`

to the covariance matrix `cov`

. This is a wrapper of `StatsBase.cor2cov.`

`AnimalBreedingTools.solve_pcg`

— Method`x = solve_pcg(A,b; maxiter=5000,eps=1e-15,reset=50,verbose=false)`

Solve the symmetric positive definite system of equations: `A`

as sparse left-hand side and `b`

as a vector of right-hand side. The Jacobi preconditioner (`diag(A)`

) is used. The maximum number of iterations is specified with `maxiter`

. Convergence criterion `eps`

is `norm(b-A*x)/norm(b)`

. The residual vector will be recalculated every `reset`

rounds. The details will be shown with `verbose=true`

.

`julia> x = solve_pcg(A,b)`

`AnimalBreedingTools.takahashi_ldlt_symtrid!`

— Method`takahashi_ldlt_symtrid!(dv,ev)`

Compute the sparse inverse (as known as selected inverse) of the LDLT factor of a symmetric trigiagonal matrix. The LDLT factor is stored in two vectors (`dv`

and `ev`

), which should be calculated with `ldlt_symtrid`

.

`AnimalBreedingTools.unvech`

— FunctionA = unvech(v; n)

Convert a half-stored vector `v`

to a square matrix `A`

with order `n`

. It is the inverse function of `vech`

. When `n`

is missing, the function will guess the order.

`AnimalBreedingTools.v2r`

— Method`corr = v2r(cov)`

Convert a covariance matrix `cov`

to the correlation matrix `corr`

. This is a wrapper of `StatsBase.cov2cor.`

`AnimalBreedingTools.vech`

— Methodv = vech(A)

Convert a square matrix to a half-stored vector. It was written by Steven G. Johnson; see `https://discourse.julialang.org/t/half-vectorization/7399`

for details.