diff --git a/blockedcholesky/.gitignore b/blockedcholesky/.gitignore new file mode 100644 index 0000000..7d2eaf7 --- /dev/null +++ b/blockedcholesky/.gitignore @@ -0,0 +1,8 @@ +/.luarc.json +/Manifest*.toml +/jds.bst +/jdsart.cls +/*.ipynb +/*.pdf +/data/ml-32m.zip +/data/ratings.arrow diff --git a/blockedcholesky/BlockedCholeskyMM.qmd b/blockedcholesky/BlockedCholeskyMM.qmd index 632f3c3..b575458 100644 --- a/blockedcholesky/BlockedCholeskyMM.qmd +++ b/blockedcholesky/BlockedCholeskyMM.qmd @@ -21,43 +21,93 @@ format: pdf-engine: pdflatex keep-tex: true runninghead: "Mixed Models in Julia" + jds-pdf: + keep-tex: true + pdf-engine: pdflatex + runninghead: "Mixed Models in Julia" + fig-width: 6 + fig-height: 3.5 + fig-format: pdf + include-in-header: + - file: texincludes.tex + jds-html: default +classoption: + - letterpaper + - review author: - - name: Douglas Bates - email: dmbates@gmail.com + - first: Douglas + first-init: D. + last: Bates + marker: '1' + footnote-id: 1 + email: bates@stat.wisc.edu orcid: 0000-0001-8316-9503 - affiliation: - - name: University of Wisconsin--Madison - city: University of Wisconsin--Madison - state: WI - url: https://www.wisc.edu - department: Department of Statistics - - name: Phillip M. Alday - email: phillip@phillipalday.com + - first: Phillip + first-init: P. + last: Alday + marker: 2 orcid: 0000-0002-9984-5745 - affiliation: - - name: Beacon Biosignals - city: Beacon Biosignals - url: https://beacon.bio - - name: Ajinkya H. Kokandakar - email: ajinkya@stat.wisc.edu - orcid: 0000-0001-6628-2272 - affiliations: - - name: University of Wisconsin--Madison - city: University of Wisconsin--Madison - state: WI - url: https://www.wisc.edu - department: Department of Statistics +# email: phillip@phillipalday.com + - first: Ajinkya + first-init: A. + last: Kokandakar + marker: 1 +footnote: + - content: Corresponding author + type: corresp + id: 1 +affiliation: + - institution: University of Wisconsin - Madison + prefix: Department of Statistics + country: U.S.A. + marker: 1 + - institution: Beacon Biosignals + country: U.S.A. + marker: 2 +# pretitle: "DISCUSSION" # optional, e.g., DISCUSSION, INVITED COMMENTARY +# received: "2024-12-01" +# revised: "2025-01-10" +# accepted: "2025-01-15" +# publishedonline: "2025-02-01" +## preamble: | +## \usepackage{...} # additional LaTeX packages if needed +# author: +# - name: Douglas Bates +# email: dmbates@gmail.com +# orcid: 0000-0001-8316-9503 +# affiliation: +# - name: University of Wisconsin--Madison +# city: University of Wisconsin--Madison +# state: WI +# url: https://www.wisc.edu +# department: Department of Statistics +# - name: Phillip M. Alday +# email: phillip@phillipalday.com +# orcid: 0000-0002-9984-5745 +# affiliation: +# - name: Beacon Biosignals +# city: Beacon Biosignals +# url: https://beacon.bio +# - name: Ajinkya H. Kokandakar +# email: ajinkya@stat.wisc.edu +# orcid: 0000-0001-6628-2272 +# affiliations: +# - name: University of Wisconsin--Madison +# city: University of Wisconsin--Madison +# state: WI +# url: https://www.wisc.edu +# department: Department of Statistics monofont: JuliaMono date: last-modified date-format: | DD MMMM YYYY [-- git #{{< git-rev short='true' >}}] abstract: | - @bates.maechler.etal:2015 described the evaluation of the profiled log-likelihood of a linear mixed-effects model by updating a sparse, symmetric positive-definite matrix and computing its Cholesky factor, as implemented in the [lme4]{.pkg} package for [R]{.proglang}. - Here we present enhancements to the derivation and theoretical presentation of the result and to its implementation using a blocked Cholesky factorization in the [MixedModels.jl]{.pkg} package for [Julia]{.proglang} [@Bezanson2017]. + @bates.maechler.etal:2015 described the evaluation of the profiled log-likelihood of a linear mixed-effects model by updating a sparse, symmetric positive-definite matrix and computing its Cholesky factor, as implemented in the {{< pkg lme4 >}} package for {{< proglang R >}}. + Here we present enhancements to the derivation and theoretical presentation of the result and to its implementation using a blocked Cholesky factorization in the {{< pkg MixedModels.jl >}} package for {{< proglang Julia >}} [@Bezanson2017]. The gain in computational efficiency is primarily due to three factors: (1) the new derivation allows us to compute the penalized residual sum of squares without computing the conditional estimates of the fixed-effects parameters and the conditional modes of the random effects at each optimization step, (2) the blocked Cholesky representation and careful ordering of the random effects terms reduces the amount of "fill-in" that occurs during the Cholesky factorization, and - (3) the multiple dispatch feature of the [Julia]{.proglang} language allows us to use specialized algorithms for different kinds of matrices instead of relying on generic algorithms during the Cholesky factorization. + (3) the multiple dispatch feature of the {{< proglang Julia >}} language allows us to use specialized algorithms for different kinds of matrices instead of relying on generic algorithms during the Cholesky factorization. To show the effectiveness of the blocked Cholesky approach we use it to fit a linear mixed model to over 32 million ratings of movies in the MovieLens `ml-32m` [@Harper2016] data set. The model incorporates random effects for over 200,000 movies and over 80,000 participants. Further enhancements to these computational methods are suggested. @@ -73,19 +123,19 @@ julia: - --project=@. --- -## Introduction {#sec-intro} +# Introduction {#sec-intro} Mixed effects models are a rich class of statistical models frequently used when the measurement units exhibit grouping structures. They are especially useful when the data exhibit non-nested grouping, for instance when there are two or more unrelated grouping variables such as measurements taken on items and users [see @baayendavidsonbates2008a]; or partially nested groupings, for instance when measuring test scores for students over time when the students can change schools [see @fuhner2021age]. -There is a long history of computational tools for fitting mixed effects models, including the `PROC MIXED` procedure in [SAS]{.proglang}, and the [nlme]{.pkg}, [lme4]{.pkg} and more recently [glmmTMB]{.pkg} packages in the [R]{.proglang} ecosystem. -The extensive use of these models is shown by the fact that @bates.maechler.etal:2015, which describes the methodology and use of [lme4]{.pkg}, is among the top 25 most-cited papers of the twenty-first century [@Pearson2025]. -We present [MixedModels.jl]{.pkg} a new implementation of mixed effects models in the [Julia]{.proglang} programming language. +There is a long history of computational tools for fitting mixed effects models, including the {{< pkg MIXED>}} procedure in {{< proglang SAS >}}, and the {{< pkg nlme >}}, {{< pkg lme4 >}} and more recently {{< pkg glmmTMB >}} packages in the {{< proglang R >}} ecosystem. +The extensive use of these models is shown by the fact that @bates.maechler.etal:2015, which describes the methodology and use of {{< pkg lme4 >}}, is among the top 25 most-cited papers of the twenty-first century [@Pearson2025]. +We present {{< pkg MixedModels.jl >}} a new implementation of mixed effects models in the {{< proglang Julia >}} programming language. This implementation relies on two innovations, one analytical and the other computational. -At the heart of [MixedModels.jl]{.pkg}, there is a new way of formulating the profiled log-likelihood evaluation that expresses this objective without requiring explicit evaluation of the conditional parameter estimates. -On the computational level, the evaluation of the objective requires the Cholesky factor of a large, sparse, positive-definite symmetric matrix, as in [lme4]{.pkg}. -This matrix has a blocked structure that is explicitly exploited in [MixedModels.jl]{.pkg} but not in [lme4]{.pkg}. +At the heart of {{< pkg MixedModels.jl >}}, there is a new way of formulating the profiled log-likelihood evaluation that expresses this objective without requiring explicit evaluation of the conditional parameter estimates. +On the computational level, the evaluation of the objective requires the Cholesky factor of a large, sparse, positive-definite symmetric matrix, as in {{< pkg lme4 >}}. +This matrix has a blocked structure that is explicitly exploited in {{< pkg MixedModels.jl >}} but not in {{< pkg lme4 >}}. More concretely, the Cholesky factor consists of dense and sparse blocks, and using specific algorithms to handle these blocks differently leads to more computationally efficient algorithms. -We exploit the multiple dispatch feature of [Julia]{.proglang} to do exactly that -- use sparse matrix algorithms for the sparse blocks, and level 3 BLAS and LAPACK subroutines for the dense blocks. +We exploit the multiple dispatch feature of {{< proglang Julia >}} to do exactly that -- use sparse matrix algorithms for the sparse blocks, and level 3 BLAS and LAPACK subroutines for the dense blocks. Together, these contributions allow us to fit linear mixed models to very large datasets. @bates.maechler.etal:2015 provided a formulation for linear mixed-effects models and showed how maximum likelihood estimates of the parameters in the model can be determined by optimizing a profiled objective function, which is negative twice the profiled log-likelihood. @@ -95,27 +145,27 @@ In this paper, we show that it is not necessary to evaluate the solution to the The objective can be evaluated from a factorization that is an intermediate step in solving the PLS problem. The model matrix, $\mbfZ$, for the random effects in such models is defined in blocks according to the "grouping factors" for random-effects terms in the model formula. -The blocked structure of $\mbfZ$ and of its cross-product matrix, $\mbfZ^\top\mbfZ$, (also called the Gram matrix of the columns of $\mbfZ$) is described in great detail in @bates.maechler.etal:2015[Section 2] when showing how a model formula and data frame create the structures used to evaluate the profiled objective. -However, the implementation of the profiled log-likelihood evaluation in the [lme4]{.pkg} package for [R]{.proglang} is not explicitly based upon the blocks --- [lme4]{.pkg} uses a more general sparse Cholesky factorization as implemented in the [CHOLMOD]{.pkg} library [@10.1145/1391989.1391995]. +The blocked structure of $\mbfZ$ and of its cross-product matrix, $\mbfZ^\top\mbfZ$, (also called the Gram matrix of the columns of $\mbfZ$) are described in great detail in @bates.maechler.etal:2015[Section 2] when showing how a model formula and data frame create the structures used to evaluate the profiled objective. +However, the implementation of the profiled log-likelihood evaluation in the {{< pkg lme4 >}} package for {{< proglang R >}} is not explicitly based upon the blocks --- {{< pkg lme4 >}} uses a more general sparse Cholesky factorization as implemented in the {{< pkg CHOLMOD >}} library [@10.1145/1391989.1391995]. -An important part in a general sparse Cholesky factorization, such as implemented in [CHOLMOD]{.pkg}, is determining a *fill-reducing permutation* of the rows and columns of the symmetric, positive-definite matrix to be factored. -The main difference between the methods in [MixedModels.jl]{.pkg} and those in [lme4]{.pkg} is that in [MixedModels.jl]{.pkg} it is the blocks that are permuted, if necessary, rather than individual rows and columns. -In other words, the particular structure of the blocks in the Gram matrix is exploited in [MixedModels.jl]{.pkg} rather than relying upon general sparse Cholesky techniques that are not able to exploit this structure. +An important part in a general sparse Cholesky factorization, such as implemented in {{< pkg CHOLMOD >}}, is determining a *fill-reducing permutation* of the rows and columns of the symmetric, positive-definite matrix to be factored. +The main difference between the methods in {{< pkg MixedModels.jl >}} and those in {{< pkg lme4 >}} is that in {{< pkg MixedModels.jl >}} it is the blocks that are permuted, if necessary, rather than individual rows and columns. +In other words, the particular structure of the blocks in the Gram matrix is exploited in {{< pkg MixedModels.jl >}} rather than relying upon general sparse Cholesky techniques that are not able to exploit this structure. The structure of the blocks in the Gram matrix depends upon the number and types of random-effects terms in the model. A general and efficient implementation of the blocked Cholesky factorization requires the ability to define operations for many different combinations of specialized forms of matrices. -[Julia]{.proglang} [@Bezanson2017] is an ideal environment for this type of calculation because of its efficient implementation of multiple dispatch, its just-in-time compilation, and its rich type system. +{{< proglang Julia >}} [@Bezanson2017] is an ideal environment for this type of calculation because of its efficient implementation of multiple dispatch, its just-in-time compilation, and its rich type system. The rest of the paper is organized as follows. In @sec-profiled we recap the derivation of the profiled log-likelihood from @bates.maechler.etal:2015, in a slightly different form from that paper, and show that the objective can be evaluated directly from the Cholesky factorization of an augmented positive-definite matrix. That is, explicit evaluation of the solution to the PLS problem is not necessary. -In @sec-structure we use an example of a fitted mixed-effects model to illustrate our implementation of the blocked Cholesky factorization in [MixedModels.jl]{.pkg}. +In @sec-structure we use an example of a fitted mixed-effects model to illustrate our implementation of the blocked Cholesky factorization in {{< pkg MixedModels.jl >}}. In @sec-application we use an example of a mixed-effects model fit to a large data set with partially crossed grouping factors for the random effects to examine how the compute time per evaluation of the objective and the memory footprint of the model depend on various dimensions of the model. We conclude with a few additional comments in @sec-further-comments followed by a short summary and discussion in @sec-summary. -## Profiled log-likelihood of linear mixed-effects models {#sec-profiled} +# Profiled log-likelihood of linear mixed-effects models {#sec-profiled} -The types of linear mixed-effects models that can be fit by [lme4]{.pkg} or [MixedModels.jl]{.pkg} are based on two vector-valued random variables: the $q$-dimensional vector of random effects, $\mcB$, and the $n$-dimensional response vector, $\mcY$. +The types of linear mixed-effects models that can be fit by {{< pkg lme4 >}} or {{< pkg MixedModels.jl >}} are based on two vector-valued random variables: the $q$-dimensional vector of random effects, $\mcB$, and the $n$-dimensional response vector, $\mcY$. The unconditional distribution of $\mcB$ and the conditional distribution of $\mcY$, given $\mcB=\mathbf{b}$, are defined as multivariate Gaussian distributions of the form $$ \begin{aligned} @@ -123,14 +173,14 @@ $$ \mcB&\sim\mcN(\mathbf{0},\mbfSigma_\theta)\ , \end{aligned} $$ {#eq-LMMdistr} - +\noindent where $\mbfX$ is the $n\times p$ model matrix for the fixed-effects terms, $\mbfbeta$ is the $p$-dimensional fixed-effects parameter vector and $\mathbf{Z}$ is the $n\times q$ model matrix for the random effects. The $q\times q$ symmetric variance-covariance matrix $\mathrm{Var}(\mcB)=\mbfSigma_\theta$ depends on the *variance-component parameter vector* $\mbftheta$ through a lower triangular *relative covariance factor* $\mbfLambda_\theta$ as $$ \mbfSigma_\theta=\sigma^2\mbfLambda_\theta\mbfLambda_\theta^\top\ . $$ {#eq-relcovfac} - +\noindent (Recall that the lower Cholesky factor is generally written $\mathbf{L}$. In this case the lower Cholesky factor contains parameters and is named with the corresponding Greek letter, $\mbfLambda$.) @@ -151,7 +201,7 @@ This occurs only with data that have been (incorrectly) simulated without error. In practice, we can safely assume that $\sigma>0$. However, the estimated $\mbfLambda_\theta$, like $\mbfSigma_\theta$, can be singular. -The computational methods in [lme4]{.pkg} and in [MixedModels.jl]{.pkg} are based on $\mbfLambda_\theta$ and do not require evaluation of $\mbfSigma_\theta$. +The computational methods in {{< pkg lme4 >}} and in {{< pkg MixedModels.jl >}} are based on $\mbfLambda_\theta$ and do not require evaluation of $\mbfSigma_\theta$. In fact, $\mbfSigma_\theta$ is explicitly evaluated only at the converged parameter estimates. To express the log likelihood without the precision matrix, we define the spherical random effects, $\mcU\sim\mcN(\mathbf{0},\sigma^2\mbfI_q)$, which in turn determine $\mcB$ as @@ -159,7 +209,7 @@ To express the log likelihood without the precision matrix, we define the spheri $$ \mcB=\mbfLambda_\theta\mcU\ . $$ {#eq-sphericalre} - +\noindent Although it may seem more intuitive to write $\mcU$ as a linear transformation of $\mcB$, we cannot do that when $\mbfLambda_\theta$ is singular, which is why @eq-sphericalre is in the form shown. We can easily verify that @eq-sphericalre provides the desired distribution for $\mcB$. @@ -171,7 +221,7 @@ $$ \mbfLambda_\theta\,\mathbb{E}\left[\mcU\right]= \mbfLambda_\theta\mathbf{0}=\mathbf{0}\ , $$ - +\noindent and covariance matrix $$ @@ -189,7 +239,7 @@ $$ \mcU&\sim\mcN(\mathbf{0},\sigma^2\mbfI_q)\ . \end{aligned} $$ {#eq-condygivenu} - +\noindent The joint density function for $\mcY$ and $\mcU$ is the product of the densities of the two distributions in @eq-condygivenu. That is @@ -256,7 +306,7 @@ r^2_\theta(\mathbf{u},\mbfbeta) \| \mathbf{r}_{Xy}-\mbfR_{XX}\mbfbeta\|^2 + r_{yy}^2\ , \end{aligned} $$ {#eq-penalized-rss} - +\noindent using the Cholesky factor of the blocked matrix, $$ @@ -281,7 +331,7 @@ $$ \end{bmatrix}\ . \end{aligned} $$ {#eq-bigCholfac} - +\noindent (In practice we evaluate the blocks of the lower Cholesky factor, such as $\mbfL_{ZZ}=\mbfR_{ZZ}^{\top}$, of this matrix but the equations are slightly easier to write in terms of the upper Cholesky factor.) Note that the block in the upper left, $\mbfLambda_\theta^{\top}\mathbf{Z^{\top}Z}\mbfLambda_\theta+\mbfI$, is positive definite even when $\mbfLambda_\theta$ is singular, because @@ -290,7 +340,7 @@ $$ \mathbf{u}^{\top}\left(\mbfLambda_\theta^{\top}\mathbf{Z^{\top}Z}\mbfLambda_\theta+\mbfI\right)\mathbf{u} = \left\|\mathbf{Z}\mbfLambda_\theta\mathbf{u}\right\|^2 +\left\|\mathbf{u}\right\|^2\ , $$ {#eq-Cholfacupperleft} -and the first term is non-negative while the second is positive if $\mathbf{u}\ne\mathbf{0}$. +\noindent and the first term is non-negative while the second is positive if $\mathbf{u}\ne\mathbf{0}$. Thus, the triangular $\mbfR_{ZZ}$, with positive diagonal elements, can be evaluated and its determinant, $\left|\mbfR_{ZZ}\right|$, which is the product of its diagonal elements, is also positive. This determinant appears in the marginal density of $\mcY$, from which the likelihood of the parameters is evaluated. @@ -300,13 +350,13 @@ To evaluate the likelihood, $$ L(\mbftheta,\mbfbeta,\sigma|\mbfy) = \int_\mathbf{u} f_{\mcY,\mcU}(\mbfy,\mathbf{u})\, d\mathbf{u}\ , $$ {#eq-likelihood-abstract} - +\noindent we isolate the part of the joint density that depends on $\mathbf{u}$ and perform a change of variable $$ \mathbf{v}=\mbfR_{ZZ}\mathbf{u}+\mbfR_{ZX}\mbfbeta-\mathbf{r}_{Zy} \ , $$ {#eq-u-system} - +\noindent with the *Jacobian matrix* $$ @@ -327,39 +377,39 @@ $$ \end{aligned} \end{aligned} $$ {#eq-likelihood-integral} - +\noindent from which we obtain the likelihood as - $$ L(\mbftheta,\mbfbeta,\sigma;\mbfy)= \frac{|\mbfR_{ZZ}|^{-1}}{(2\pi\sigma^2)^{n/2}} \exp\left(-\frac{r_{yy}^2 + \|\mbfR_{XX}(\mbfbeta-\widehat{\mbfbeta})\|^2}{2\sigma^2}\right)\ , $$ {#eq-likelihood} - +\noindent where the conditional estimate, $\widehat{\mbfbeta}$, given $\mbftheta$, satisfies $$ \mbfR_{XX}\widehat{\mbfbeta} = \mathbf{r}_{Xy}\ . $$ - +\noindent Setting $\mbfbeta=\widehat{\mbfbeta}$ and taking the logarithm provides the estimate of $\sigma^2$, given $\mbftheta$, as - $$ -\widehat{\sigma^2}=\frac{r_yy^2}{n}\ , +\widehat{\sigma^2}=\frac{r_{yy}^2}{n}\ , $$ {#eq-sigma-hat} - +\noindent which gives the *profiled log-likelihood*, $\ell(\mbftheta|\mbfy)=\log L(\mbftheta,\widehat{\mbfbeta},\widehat{\sigma})$, on the deviance scale, as $$ -2\ell(\mbftheta|\mbfy)=2\log(|\mbfR_{ZZ}|) + -n\left(1+\log\left(\frac{2\pi r_{yy}^2(\mbftheta)}{n}\right)\right)\ . +n\left(1+\log\left(\frac{2\pi r_{yy}^2(\mbftheta)}{n}\right)\right)\ , $$ {#eq-profiled-log-likelihood} +\noindent +which is an affine function of the logarithms of the diagonal elements of $\mbfR$. One of the interesting aspects of this formulation is that it is not necessary to solve for the conditional estimate of $\mbfbeta$ or the conditional modes of the random effects when evaluating the log-likelihood. -The two values needed for the log-likelihood evaluation, $2\log(|\mbfR_{ZZ}|)$ and $r_\mathbf{yy}^2$, are obtained directly from the diagonal elements of the Cholesky factor. +The two values needed for the log-likelihood evaluation, $2\log(|\mbfR_{ZZ}|)$ and $r_{yy}^2$, are obtained directly from the diagonal elements of the Cholesky factor. Furthermore, $\mbfOmega_{\theta}$ and, from that, the Cholesky factor, $\mbfR_{\theta}$, and the objective to be optimized can be evaluated for a given value of $\mbftheta$ from the Gram matrix of the columns of $\mbfZ$, $\mbfX$ and $\mbfy$, @@ -370,34 +420,51 @@ $$ \mbfy^{\top}\mathbf{Z} & \mbfy^{\top}\mbfX & \mbfy^{\top}\mbfy \end{bmatrix}\ , $$ {#eq-A} - +\noindent and $\mbfLambda_{\theta}$. -In [MixedModels.jl]{.pkg} the `LinearMixedModel` struct contains a blocked representation of $\mathbf{A}$ in the `A` field and a similarly structured lower-triangular blocked array in the `L` field. -Evaluation of the objective simply involves updating a representation of the relative covariance factor $\mbfLambda_\theta$ from $\mbftheta$, forming $\mbfOmega_\theta$ (@eq-bigCholfac) from $\mbfA$ and $\mbfLambda_\theta$ then evaluating its lower Cholesky factor $\mbfL_\theta$. +In {{< pkg MixedModels.jl >}} the `LinearMixedModel` struct contains a blocked representation of $\mathbf{A}$ in the `A` field and a similarly structured lower-triangular blocked array in the `L` field. +Evaluation of the objective simply involves updating a representation of the relative covariance factor $\mbfLambda_\theta$ from $\mbftheta$, forming $\mbfOmega_\theta$ (@eq-bigCholfac) from $\mbfA$ and $\mbfLambda_\theta$, then evaluating its lower Cholesky factor, $\mbfL_\theta$. In @eq-A the blocking structure is according to all the random effects, followed by the fixed-effects followed by the response. However, as described in @bates.maechler.etal:2015[Section 3], the structure of $\mbfZ$, $\mbfZ^{\top}\mbfZ$ and $\mbfLambda_\theta$ is further subdivided into blocks according to the random-effects terms in the model formula, which we cover in the next section. -## Random-effects terms in mixed-model formulas {#sec-structure} +# Random-effects terms in mixed-model formulas {#sec-structure} -The user interface in [MixedModels.jl]{.pkg} is intentionally similar to that of [lme4]{.pkg}. +The user interface in {{< pkg MixedModels.jl >}} is intentionally similar to that of {{< pkg lme4 >}}. A model is specified, as described in @bates.maechler.etal:2015[section 2.1], by a formula of the form ```julia @formula(resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...) ``` +\noindent "where `FEexpr` is an expression determining the columns of the fixed-effects model matrix, $\mbfX$, and the random-effects terms, `(REexpr1 | factor1)`, `(REexpr2 | factor2)`, etc. together determine both the random-effects model matrix $\mbfZ$ and the structure of the relative covariance factor $\mbfLambda_\theta$." -(In [Julia]{.proglang} a formula expression such as this must be enclosed in a call to the `@formula` macro.) +(In {{< proglang Julia >}} a formula expression such as this must be enclosed in a call to the `@formula` macro.) The expression on the right hand side of the vertical bar, `|`, in a random-effects term is called the *grouping factor* for those random effects. -### From formula to model matrices +## From formula to model matrices The process of converting a mixed-model formula for a data set into the random-effects model matrix, $\mbfZ$, the Gram matrix, $\mbfZ^{\top}\mbfZ$, and the mapping from the relative covariance parameter, $\mbftheta$, to the relative covariance factor, $\mbfLambda_\theta$, is described in considerable detail in @bates.maechler.etal:2015[Sections 2.2 and 2.3]. In particular, Tables 3 and 4 of @bates.maechler.etal:2015 provide the dimensions and the names of some of the components of these matrices and should be kept close at hand when reading this section. -For convenience, Table 3 from @bates.maechler.etal:2015 has been reproduced with slight changes as @tbl-re-sizes in Appendix -@sec-re-dim. +For convenience, Table 3 from @bates.maechler.etal:2015 has been reproduced with slight changes as @tbl-re-sizes. + + +| Symbol | Description | `insteval` example | +|:--------------------|:--------------------------------------------------------------|:---------------------------| +| $n$ | length of the response vector $\mcY$ | $n = 73{,}421$ | +| $p$ | no. of fixed effects in the model matrix $\mbfX$ | $p = 2$ | +| $k$ | no. of random effects terms | $k = 3$ (`s`, `d`, `dept`) | +| $p_i$ | no. of fixed effects terms in the $i^\text{th}$ random effect | $(1,1,2)$ | +| $\ell_i$ | no. of levels of the $i^\text{th}$ grouping factor | $(2972, 1128, 14)$ | +| $q_i = p_i \ell_i$ | no. of columns of the term-wise model matrix | $(2972, 1128, 28)$ | +| $q = \sum_{i}q_i$ | no. of columns in the random effects model matrix $\mbfZ$ | $q = 4{,}128$ | +| $q + p + 1$ | no. of columns (and rows) of the augmented matrix $\mbfA$ | $4131$ | + + +: Sizes of model components {#tbl-re-sizes tbl-colwidths="[18,62,20]"} + To recap some of the material in @bates.maechler.etal:2015[Sections 2.2 and 2.3], in a model with $k$ random-effects terms, an $n\times p_i$ model matrix $\mbfX_i, i=1,\dots,k$ is evaluated from the expression `REterm`$_i$ and the data. The $\mbfX_i, i=1,\dots,k$ are (implicitly) combined with the indicator matrices, $\mathbf{J}_i, i=1,\dots,k$ of the $i$th grouping factor to produce the sparse model matrix blocks $\mbfZ_i, i=1,\dots,k$. @@ -407,44 +474,37 @@ Generally the $p_i, i=1,\dots,k$ are small. Random effects terms of the form `(1|g)`, i.e. $p_i=1$, are not uncommon. In such cases, called simple, scalar random-effects terms, the matrix $\mbfZ_i$ is exactly the indicator matrix $\mathbf{J}_i$ and its $\ell_i\times\ell_i$ Gram matrix, $\mbfZ_i^{\top}\mbfZ_i$, is diagonal with the incidence counts of the levels on the diagonal. -The blocked derivation in @bates.maechler.etal:2015 carries over into the model representation in [MixedModels.jl]{.pkg} with one modification --- in [lme4]{.pkg} the (implicit) partitioning of $\mbfZ$ into $k$ vertical blocks is according to the random-effects terms in the formula, whereas in [MixedModels.jl]{.pkg} it is according to the grouping factors. +The blocked derivation in @bates.maechler.etal:2015 carries over into the model representation in {{< pkg MixedModels.jl >}} with one modification --- in {{< pkg lme4 >}} the (implicit) partitioning of $\mbfZ$ into $k$ vertical blocks is according to the random-effects terms in the formula, whereas in {{< pkg MixedModels.jl >}} it is according to the grouping factors. If more than one random-effects terms have the same grouping factor they are amalgamated in the model representation. Making the number of blocks small reduces the complexity of the evaluation of the profiled likelihood while retaining simplicity in the factorization of the first column of blocks (see [@sec-one-one-block]). ### An example using instructor evaluation data -The `InstEval` dataset from the [lme4]{.pkg} package is reproduced as `insteval` in [MixedModels.jl]{.pkg}. -(The [Julia]{.proglang} convention is to use lower-case for names of objects and camel-case for names of types and for package names.) +The `InstEval` dataset from the {{< pkg lme4 >}} package is reproduced as `insteval` in {{< pkg MixedModels.jl >}}. +(The {{< proglang Julia >}} convention is to use lower-case for names of objects and camel-case for names of types and for package names.) It consists of nearly 75,000 evaluations (`y`) of instructors (`d` for *Dozent*, German for "instructor") at ETH-Zürich by students (`s`). Other covariates, such as the department (`dept`) that offered the course and whether or not it was a `service` course, are included. -Further details are given on the [R]{.proglang} help page for `lme4::InstEval`. +Further details are given on the {{< proglang R >}} help page for `lme4::InstEval`. -In [Julia]{.proglang} we attach the packages to be used, load the data set, convert it to a `DataFrame` and evaluate some descriptive statistics on its columns to produce @tbl-insteval-summary, +In {{< proglang Julia >}} we attach the packages to be used, load the data set, convert it to a `DataFrame` and evaluate some descriptive statistics on its columns to produce @tbl-insteval-summary, which shows that a total of 2972 students evaluated a total of 1128 instructors in courses from 14 different departments. About 43% of the evaluations were from service courses. ```{julia} -#| output: false -using AlgebraOfGraphics, CairoMakie, PrettyTables, Printf, SparseArrays, LinearAlgebra -using Arrow, Chairmarks, CSV, DataFrames, MixedModels -dat = DataFrame(MixedModels.dataset("insteval")) -dat.service = float.(dat.service .== "Y") # convert to 0/1 encoding -sumry = describe(dat, :min, :max, :mean, :nunique, :eltype) -``` - -```{julia} -#| echo: false -#| output: asis -#| caption: summary #| tbl-cap: Summary of the `insteval` dataset #| label: tbl-insteval-summary -fmt_nothing = (v, i, j) -> isnothing(v) ? "-" : v -pretty_table( - sumry, - formatters = (ft_round(2, [4]), fmt_nothing), - show_subheader = false, backend = Val(:latex)) +using AlgebraOfGraphics, CairoMakie, Printf, SparseArrays +using Arrow, Chairmarks, CSV, DataFrames, LinearAlgebra, MixedModels +dat = DataFrame(MixedModels.dataset("insteval")) +dat.service = (dat.service .== "Y") # convert to Boolean values +show( + describe(dat, :min, :max, :mean, :nunique, :eltype); + row_labels=nothing, + eltypes=false, + summary=false, +) ``` ```{julia} @@ -470,8 +530,8 @@ The model formula is chosen to illustrate certain aspects of the conversion from In particular, we choose independent random effects for `dept` and for `service` by `dept` to show that distinct terms with the same grouping factor will be amalgamated in the internal representation. ```{julia} -form = @formula(y ~ 1 + service + (1|d) + (1|s) + (1|dept) + (0 + service|dept)) -m1 = fit(MixedModel, form, dat, progress=false) # suppress display of a progress bar +form = @formula(y ~ 1 + service + (1|d) + (1|s) + (1|dept) + (0+service|dept)) +m1 = fit(MixedModel, form, dat, progress=false) # suppress display of progress bar print(m1) ``` @@ -479,20 +539,20 @@ Notice that, in the formula, the random effects term for instructor, `(1|d)`, oc In the model the `s` term occurs first because the number of random effects from that term, 2972, is greater than that from `d`, 1128. Also the two terms for `dept` are amalgamated into a single term in the model. -The optimization of the profiled objective required about 175 evaluations of the objective function and took about a second on a laptop computer, as shown in this benchmark using the `@be` macro from the [Chairmarks.jl]{.pkg} package. +The optimization of the profiled objective required about 175 evaluations of the objective function and took about a second on a laptop computer, as shown in this benchmark using the `@be` macro from the {{< pkg Chairmarks.jl >}} package. ```{julia} @be fit($MixedModel, $form, $dat; progress=false) seconds=8 ``` -The environment in which these benchmarks were performed is described in @sec-app-env. +The environment in which these benchmarks were performed is described in the supplemental materials. Note that we use the Apple Accelerate implementation of the BLAS/Lapack routines on computers with Apple M-series processors and Intel's MKL BLAS implementation on computers with x86_64 processors. -In Julia hot-swapping BLAS implementations is as easy as loading the appropriate package, such as [AppleAccelerate.jl]{.pkg} or [MKL.jl]{.pkg}. +In Julia hot-swapping BLAS implementations is as easy as loading the appropriate package, such as {{< pkg AppleAccelerate.jl >}} or {{< pkg MKL.jl >}}. Individual evaluations of the profiled objective, including installing a new covariance parameter vector, updating the blocked Cholesky factorization, and evaluating @eq-profiled-log-likelihood, takes a few milliseconds on the same laptop computer. ```{julia} -@be objective(updateL!(setθ!($m1, $(m1.θ)))) +@be objective!($m1, $(m1.optsum.final)) ``` The random effects terms define scalar random effects for `s` and `d` and vector-valued random effects for `dept` in this model, as reflected in the data type of each of the `reterms` of the model @@ -500,60 +560,60 @@ The random effects terms define scalar random effects for `s` and `d` and vector ```{julia} show([typeof(rt) for rt in m1.reterms]) ``` - -In the notation of @bates.maechler.etal:2015[Section 2] the second parameter of each of these types is $p_i, i = 1,\dots,3$. +\noindent +In the notation of @bates.maechler.etal:2015[Section 2] the second parameter of each of these types is $p_i, i = 1,2,3$. The number of levels of the grouping factors are ```{julia} show([length(rt.levels) for rt in m1.reterms]) ``` +\noindent +In the notation of @bates.maechler.etal:2015[Section 2] these are $\ell_i, i=1,2,3$. -In the notation of @bates.maechler.etal:2015[Section 2] these are $\ell_i, i=1,\dots,3$. - -Finally the number of columns in each of the $\mbfZ_i, i=1,\dots,3$ is +Finally the numbers of columns in each of the $\mbfZ_i, i=1,2,3$ are ```{julia} show([size(rt, 2) for rt in m1.reterms]) ``` - +\noindent For a summary of the model component sizes see @tbl-re-sizes. -The blocks of the augmented Gram matrix, $\mbfA$, and the lower Cholesky factor, $\mbfL_\theta=\mbfR^{\top}_\theta$, are summarized as +The blocks of the augmented Gram matrix, $\mbfA$, and the lower Cholesky factor, $\mbfL_\theta=\mbfR^{\top}_\theta$, are summarized in @tbl-blk_descrip ```{julia} -#| tbl-cap: Block Structure of the `A` matrix \label{tab:blk-descrip} +#| tbl-cap: Block Structure of the `A` and `L` matrices in the `insteval` example +#| label: tbl-blk_descrip #| tbl-cap-location: top BlockDescription(m1) ``` - +\noindent (Only the lower triangle of the symmetric matrix, $\mbfA$, is stored.) In the theoretical development the blocking structure of the Gram matrix, $\mbfA$, and the Cholesky factor (@eq-bigCholfac) is according to random effects, fixed effects and response. In practice the response vector, $\mbfy$, is concatenated (`hcat`) to the fixed effects model matrix, $\mbfX$, and the random effects model matrix, $\mbfZ$, is partitioned according to the grouping factors. -In the block structure summary if a block of $\mbfL$ and its corresponding block of $\mbfA$ have the same structure, the name of that structure is given once. +In the block structure summary (@tbl-blk_descrip) if a block of $\mbfL$ and its corresponding block of $\mbfA$ have the same structure, the name of that structure is given once. This is the case for all of the blocks except the [2,2] block, which is diagonal in $\mbfA$ but dense in $\mbfL$, the [2,3] block, which is sparse in $\mbfA$ but dense in $\mbfL$, and the [3,3] block, which is block-diagonal in $\mbfA$ and dense in $\mbfL$. As described in @bates.maechler.etal:2015[Section 3] the relative covariance factor, $\mbfLambda_\theta$ is block diagonal with the same blocking structure as the random effects parts of $\mbfA$ or $\mbfL$. The $i$th diagonal block, $i=1,\dots,k$, consists of $q_i$ repetitions on the diagonal of a $\ell_i\times\ell_i$ lower-triangular template matrix $\mbfT_i$. For a scalar random effects term (i.e. $\ell_i=1$) its block in $\mbfLambda_\theta$ is a multiple of the $q_i\times q_i$ identity matrix. -The blocks of $\mbfLambda_\theta$ are never instantiated in the [MixedModels.jl]{.pkg} implementation. There is sufficient flexibility in the multiple dispatch model in [Julia]{.proglang} that pre-multiplication by $\mbfLambda_\theta^\top$ or post-multiplication by $\mbfLambda_\theta$ can be performed from the template matrices $\mbfT_i, i=1,\dots,k$ and the `refs` property of the grouping factor alone. +The blocks of $\mbfLambda_\theta$ are never instantiated in the {{< pkg MixedModels.jl >}} implementation. There is sufficient flexibility in the multiple dispatch model in {{< proglang Julia >}} that pre-multiplication by $\mbfLambda_\theta^\top$ or post-multiplication by $\mbfLambda_\theta$ can be performed from the template matrices $\mbfT_i, i=1,\dots,k$ and the `refs` property of the grouping factor alone. The value of $\hat{\mbftheta}$ for model `m1` is ```{julia} -θ̂ = m1.θ -show(θ̂) +show(m1.optsum.final) ``` - +\noindent which maps to the template matrices ```{julia} -m1.λ +m1.lambda ``` -### Special role of the [1,1] block {#sec-one-one-block} +## Special role of the [1,1] block {#sec-one-one-block} The reason for ordering the random-effects terms so that the [1,1] block of $\mbfA$ is as large as possible is because the diagonal or block diagonal structure of this block is preserved in the [1,1] block of $\mbfL$. To illustrate the importance of choosing the first block of the augmented matrix $\mbfA$, we fit the same model with (a) the [1,1] block corresponding to the random effect term for the instructor (`1 | d`) with 1128 levels, and (b) the [1,1] block corresponding to the random effect term for students (`(1 | s)` with 2972 levels). @@ -578,7 +638,8 @@ copyto!(m2.L, L); # set L ```{julia} #| echo: false -refit!(m2; progress = false); +#| output: false +refit!(m2; progress = false) ``` ```{julia} @@ -600,7 +661,7 @@ f #| fig-cap: Sparsity pattern of the lower triangular Cholesky factor when the order of the random effects is not carefully chosen (left) and when the order is carefully chosen (right). #| label: fig-blk -second(x) = x[2] +second(x) = x[2] # why not just use `last`? reszfast = second.(size.(m1.reterms)) reszslow = second.(size.(m2.reterms)) szl, _ = size(sparseL(m1; full = true)) @@ -689,20 +750,20 @@ CairoMakie.activate!(; type="pdf"); ![Sparsity pattern of the lower triangular Cholesky factor when the order of the random effects is not carefully chosen (left) and when the order is carefully chosen (right).](lmm_side_by_s.png){#fig-blk width=100% fig-pos='!h'} --> -@fig-a shows the non-zero elements of the augmented matrix, $\mbfA$, when (a) the first block corresponsds to the random effect terms for the instructor (`(1 | d)` with 1128 levels), and (b) the first block corresponds to the student (`(1 | s)` with 2972 levels), and @fig-blk shows the non-zero elements for the corresponding lower-triangular Cholesky factor. -When performing a block-wise Cholesky factorization, the diagonal structure of the $[1,1]$ block from the input matrix is preserved in the Cholesky factor, but the diagonal structure for the subsequent diagonal blocks (e.g [2,2]) is not preserved. +@fig-a shows the non-zero elements of the augmented matrix, $\mbfA$, when (a) the first block corresponds to the random effect terms for the instructor (`(1 | d)` with 1128 levels), and (b) the first block corresponds to the student (`(1 | s)` with 2972 levels), and @fig-blk shows the non-zero elements for the corresponding lower-triangular Cholesky factor. +When performing a block-wise Cholesky factorization, the diagonal structure of the $[1,1]$ block from the input matrix is preserved in the Cholesky factor, but the diagonal structure for the subsequent diagonal blocks (e.g. [2,2]) is not preserved +(except in the case that the levels of the first grouping factor are strictly nested within the levels of the second grouping factor). While the models are theoretically identical, the latter is computationally more efficient as the largest diagonal block preserves its structure. -In fact, for the `insteval` dataset, the former order leads to approximately 4.5 million non-zero entries in the Cholesky factor, whereas choosing the random effect with the largest number of levels leads to a Cholesky factor with approximately 775,000 non-zero elements, which is an order of magnitude lower. -This translates to an order of magnitude difference between the runtimes for fitting the two versions of the model. -We demonstrate this with the `insteval` dataset in Appendix -@sec-app-re where fitting the model with the two blocks swapped is almost twenty times slower. +In fact, for the `insteval` dataset, the former order leads to approximately 4.5 million non-zero entries in the Cholesky factor, whereas choosing as first the random effect with the largest number of levels leads to a Cholesky factor with approximately 775,000 non-zero elements, which is an order of magnitude lower. +In practice this change gives an order of magnitude difference between the runtimes for fitting the two versions of the model for the `insteval` data. Because the [1,1] block of $\mbfL$ is diagonal and the [1,1] block of $\mbfLambda_\theta$ is a multiple of the identity, the remaining blocks in the first column of blocks of $\mbfL$ will have exactly the same sparsity structure as the corresponding blocks of $\mbfA$. For example, the [2,1] block of $\mbfA$ is stored as a sparse matrix because only about 2% (73421) of the $2792\times1128=3352416$ elements of that matrix are non-zero. The [2,1] block of $\mbfL$ has exactly the same sparsity structure because each element of that block of $\mbfL$ is a scalar multiple of the corresponding element of that block of $\mbfA$. -## Application to data from an observational study {#sec-application} +# Application to data from an observational study {#sec-application} -We consider an observational study -- ratings of movies by users of [movielens.org](https://movielens.org) made available at the [grouplens.org download site](https://grouplens.org/datasets/movielens) [@Harper2016] as *MovieLens 32M Dataset* (`ml-32m`), consisting of roughly 32 million ratings of over 80,000 movies by over 200,000 users. +We consider a large observational study -- ratings of movies by users of [movielens.org](https://movielens.org) made available at the [grouplens.org download site](https://grouplens.org/datasets/movielens) [@Harper2016] as the *MovieLens 32M Dataset* (`ml-32m`), consisting of roughly 32 million ratings of over 80,000 movies by over 200,000 users. The purpose of this section is to study which dimensions of the data have the greatest effect on the amount of memory used to represent the model and the time required to fit a model. To do so we use a simple model: @@ -710,16 +771,18 @@ To do so we use a simple model: frm = @formula(rating ~ 1 + (1 | userId) + (1 | movieId)) fit(LinearMixedModel, frm, dataset) ``` +\noindent where `dataset` will be different subsets of the full data explained below. -### Structure of the data +## Structure of the data The `ratings.csv` file in the `ml-32m.zip` file downloaded from the Grouplens site provides a table of about 32 million ratings with the corresponding movie number and user number. Creating a `ratings` data frame from this file and counting the number of ratings for each movie id provides the `movies` data frame with columns `movieId` and `mnrtngs`. Similarly counting the number of ratings for each user provides the `users` data frame with columns `userId` and `unrtngs`. Joining these two derived data frames with the original `ratings` data frame (on `movieId` and `userId`, respectively) allows for the `ratings` data frame to be subsetted according to a minimum number of ratings per user and/or per movie. +The supplemental materials contain a <<{ proglang Julia }>> script, `ml32-download.jl` to perform the downloading and extraction of the csv file and creation of the subsequent tables. ```{julia} #| echo: false @@ -733,7 +796,7 @@ The data from this observational study are extremely unbalanced with respect to ```{julia} #| echo: false -#| fig-cap: "Empirical distribution plots of the number of ratings per movie and per user. The horizontal axes are on a logarithmic scale." +#| fig-cap: "Empirical cumulative distribution plots of the number of ratings per movie and per user. The horizontal axes are on a logarithmic scale." #| label: fig-nrtngsecdf #| warning: false let @@ -744,7 +807,7 @@ let yminorticksvisible = true xminorticks = IntervalsBetween(10) ylabel = "Relative cumulative frequency" - nrtngs = sort(movies.nrtngs) + nrtngs = sort(movies.mnrtngs) ecdfplot( f[1, 1], nrtngs; @@ -759,7 +822,7 @@ let yminorticksvisible, ), ) - urtngs = sort(users.nrtngs) + urtngs = sort(users.unrtngs) ecdfplot( f[1, 2], urtngs; @@ -796,24 +859,46 @@ One way of dealing with the extreme imbalance in the number of observations per ### Models fit with lower bounds on ratings per user and per movie {#sec-lrgobsmods} We fit a simple model to this dataset using different thresholds on the number of ratings per movie and the number of ratings per user. -These fits were performed on a desktop computer with a generous amount of memory (96 GiB) and using 8 threads and the Intel Math Kernel Library (MKL) for numerical linear algebra. -The results are summarized in @tbl-sizespeed in Appendix -@sec-app-ml. +These fits were performed on a cloud instance with a generous amount of memory (128 GiB) and using 12 threads and the Intel Math Kernel Library (MKL) for numerical linear algebra. +The results are summarized in @tbl-sizespeed. ```{julia} #| echo: false -#| output: false sizespeed = CSV.read("./data/sizespeed.csv", DataFrame; downcast=true); ``` -#### Dimensions of the model versus cut-off values -As shown in the first panel of @fig-nratingsbycutoff, the number of ratings varies from a little over 21 million to 32 million. +```{=latex} +\ifXeTeX +\footnotesize +\else +\small +\fi +``` + +```{julia} +#| echo: false +#| tbl-cap: Summary of results for the `ml-32m` data modeling run speeds. The `mc` column is the movie cutoff, `uc` is the user cutoff, `nratings`, `nusers` and `nmvie` are, respectively, the number of observations, users and movies in the reduced data set, `modelsz` and `L22sz` are the model size and `L[2,2]` size (in GiB), `nv` is the number of objective function evaluations until convergence, `fittime` and `mnevtm` are the number of seconds required to fit the model and the mean time to perform a single objective function evaluation. +#| label: tbl-sizespeed +sizespeed = CSV.read("./data/sizespeed.csv", DataFrame; downcast=true) +show(sizespeed; row_labels=nothing, eltypes=false, allrows=true, summary=false) +``` + +```{=latex} +\normalsize +\clearpage +``` + +### Dimensions of the model versus cut-off values + +As shown in the first panel of @fig-nratingsbycutoff, the number of ratings varies from a little over 27 million to 32 million. For this set of cutoff values, the user cutoff has more impact on the number of ratings in the reduced dataset than does the movie cutoff. +A glance at @tbl-sizespeed shows that the number of users is essentially a function of only the user cutoff. ```{julia} #| echo: false #| label: fig-nratingsbycutoff -#| fig-cap: "Number of ratings (left) and movies (right) in reduced table by movie cutoff and by user cutoff" +#| fig-cap: "Number of ratings (left) and movies (right) in reduced table by movie cutoff and by user cutoff. There are three lines in the right panel but they are essentially superimposed." #| warning: false function siformat(x) @@ -845,13 +930,13 @@ draw( ) ``` -A glance at @tbl-sizespeed shows that the number of users is essentially a function of only the user cutoff. - The second panel in @fig-nratingsbycutoff shows the similarly unsurprising result that the number of movies in the reduced table is essentially determined by the movie cutoff. +(There are three different, but nearly coincident, lines in this panel.) -#### Memory footprint of the model representation {#sec-lrgobsmemprint} +### Memory footprint of the model representation {#sec-lrgobsmemprint} -@fig-memoryfootprint shows that the overall memory size (the "memory footprint") of the model representation depends primarily on the movie cutoff. (There are three lines on @fig-memoryfootprint, corresponding to the three user cutoff values, but they are nearly coincident.) +@fig-memoryfootprint shows that the overall memory size (the "memory footprint") of the model representation depends primarily on the movie cutoff. +(There are also three distinct, but nearly coincident, lines in @fig-memoryfootprint, corresponding to the three user cutoff values.) ```{julia} #| echo: false @@ -875,7 +960,7 @@ Even when requiring a minimum of 80 ratings per movie, `L[2,2]` accounts for abo ```{julia} #| echo: false -#| fig-cap: Proportion of memory footprint of the model in L[2,2] versus the overall model size (GiB). +#| fig-cap: Proportion of memory footprint of the model in L[2,2] versus the overall model size (GiB). The dots are at different movie cutoff (mc) values, with the largest memory footprints at mc=1 and the smallest at mc=50. #| label: fig-l22prop #| warning: false draw( @@ -892,9 +977,9 @@ draw( The fact that the `[2,2]` block of `A` is diagonal but the `[2,2]` block of `L` is dense lower-triangular is described as "fill-in" and leads to a somewhat unintuitive conclusion. The memory footprint of the model representation depends strongly on the number of movies, less strongly on the number of users and almost not at all on the number of ratings ([@fig-memoryfootprint]). -Note that although the dimensions of the `L[2,1]` block are larger than those of the `L[2,2]` block its memory footprint is smaller than that of the `[2,2]` block, because it is stored as a sparse matrix. +Note that although the dimensions of the `L[2,1]` block are larger than those of the `L[2,2]` block its memory footprint is smaller than that of the `L[2,2]` block, because it is stored as a sparse matrix. The matrix is about 98% zeros or, equivalently, a little over 2% nonzeros, which makes the sparse representation much smaller than the dense representation. -Also, the number and positions of the non-zeros in `L[2,1]` is exactly the same as those of in `A[2,1]`. +Also, the number and positions of the non-zeros in `L[2,1]` is exactly the same as those in `A[2,1]`. That is, there is no fill-in in the `[2,1]` block In a sense this is good news because the amount of storage required for the `[2,2]` block can be nearly cut in half by taking advantage of the fact that it is a triangular matrix. @@ -902,41 +987,41 @@ See @sec-RFP for further discussion on one approach using the *rectangular full In general, for models with scalar random effects for two incompletely crossed grouping factors, the memory footprint depends strongly on the smaller of the number of levels of the grouping factors, less strongly on the larger number of levels, and almost not at all on the number of observations. -#### Speed of log-likelihood evaluation +### Speed of log-likelihood evaluation The time required to fit a model to large data sets is dominated by the time required to evaluate the log-likelihood during the optimization of the parameter estimates. -The time for one evaluation is given in the `time per eval` column of @tbl-sizespeed. -Also given is the number of evaluations to convergence, `n eval`, and the total time to fit the model, `time`. -The reason for considering `evtime` in addition to `fittime` and `n eval` is because the `time per eval` for one model, relative to other models, is reasonably stable across computers whereas `n eval`, and hence, `time`, can be affected by seemingly trivial variations in function values resulting from different implementations of low-level calculations, such as the BLAS (Basic Linear Algebra Subroutines). +The mean time per evaluation is given in the `mnevtm` column of @tbl-sizespeed. +Also given is the number of evaluations to convergence, `nv`, and the total time to fit the model, `fittime`. +The reason for considering `mnevtm` in addition to `fittime` and `nv` is because the `mnevtm` for one model, relative to other models, is reasonably stable across computers whereas `nv`, and hence, `fittime`, can be affected by seemingly trivial variations in function values resulting from different implementations of low-level calculations, such as the BLAS (Basic Linear Algebra Subroutines). (Floating point arithmetic is not guaranteed to be associative and custom BLAS implementations may rearrange the order of operations when evaluating products and sums.) -That is, we can't expect to reproduce `n eval` exactly when fitting the same model on different computers or with slightly different versions of software; but we can expect the pattern in `time per eval` with respect to the user and movie cutoffs to be reproducible. - -As shown in @fig-evtimevsl22 the evaluation time for the objective is predominantly a function of the size of the `[2,2]` block of `L`. +That is, we can't expect to reproduce `nv` exactly when fitting the same model to the same data but on different computers or with slightly different versions of software. +However, we can expect the pattern in `mnevtm` with respect to the user and movie cutoffs to be reproducible. +As shown in @fig-mnevtmvsl22 the evaluation time for the objective is predominantly a function of the size of the `[2,2]` block of `L`. ```{julia} #| echo: false #| fig-cap: "Evaluation time for the objective (s) versus size of the [2,2] block of L (GiB)" -#| label: fig-evtimevsl22 +#| label: fig-mnevtmvsl22 #| warning: false draw( data(sizespeed) * mapping( :L22sz => "Size of [2,2] block of L (GiB)", - :evtime => "Time for one evaluation of objective (s)", + :mnevtm => "Time for one evaluation of objective (s)", color = :uc => renamer([20 => "≥20 ratings", 40 => "≥40 ratings", 80 => "≥80 ratings"]) => "per user",) * visual(ScatterLines); legend = (; position = :bottom, titleposition = :left), ) ``` -#### Reducing the amount of storage required {#sec-RFP} +### Reducing the amount of storage required {#sec-RFP} As shown in @fig-l22prop the `L[2,2]` block, when stored as a dense square matrix, can be a major contributor to the amount of storage required for the model. Because only the lower triangle of this dense square matrix is used, about half of this storage is never accessed. The rectangular full-packed (RFP) storage format [@lawn199] allows for a triangular matrix to be packed into (roughly) half the amount of storage required for the dense square matrix. -The [RectangularFullPacked.jl]{.pkg} package allows access to this format in [Julia]{.proglang}. +The {{< pkg RectangularFullPacked.jl >}} package allows access to this format in {{< proglang Julia >}}. The savings in storage comes at the expense of a more complicated algorithm for mapping matrix coordinates to storage locations. @@ -949,10 +1034,10 @@ In practical examples, $\mbfL[2,1]$ is sparse and the update is performed by ite Each pair of nonzeros in each column of $\mbfL[2,1]$ generates a (scalar) update of one of the elements of $\mbfL[2,2]$ and this is where the more expensive indexing system in the RFP formulation (relative to storing a square dense matrix) takes its toll. Our current implementation using RFP takes three to four times as long to perform the rank-k update than to do the factorization to obtain $\mbfL[2,2]$. -#### A[2,1] as a biadjacency matrix {#sec-biadjacency} +### A[2,1] as a biadjacency matrix {#sec-biadjacency} In both of the examples we have considered, the first and second random-effects terms are simple scalar terms and each observation corresponds to a unique combination of levels of these two factors. -For example, in the movielens data each rating corresponds to a unique combination of a user and a movie. +For example, in the `movielens` data each rating corresponds to a unique combination of a user and a movie. In cases like these, the $\mbfA[2,1]$ matrix is the biadjacency matrix for a bipartite graph $(U,V,E)$ where the nodes $U$ are the individual users, the nodes $V$ are the individual movies and the edges, $E$, link a movie to a user who rated it. The total number of edges in the bipartite graph is $n$, the number of observations in the data set. @@ -960,7 +1045,7 @@ Furthermore, the diagonal elements in the diagonal matrices $\mbfA[1,1]$ and $\m They are also the column- and row-sums of $\mbfA[2,1]$. It may be possible to derive new computational approaches from these relationships using graph theory, an avenue we hope to explore in future research. -## Further Comments {#sec-further-comments} +# Further Comments {#sec-further-comments} Throughout the discussion we have considered maximum likelihood estimation, which is the default estimation criterion. However, it is also possible to fit mixed models using the restricted maximum likelihood (REML) criterion by setting the optional argument `REML=true` in the call to the `fit` function. @@ -971,200 +1056,61 @@ For instance: ```{julia} mreml = fit( MixedModel, - @formula(y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)), + @formula(y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service|dept)), dat, REML = true, progress=false ) print(mreml) ``` +\noindent +The estimates obtained above are the same as those obtained from `lme4::lmer` function in {{< proglang R >}}. -The estimates obtained above are the same as those obtained from `lme4::lmer` function in [R]{.proglang}. - -We used the formula `y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)` to demonstrate that [MixedModels.jl]{.pkg} amalgamates random effects terms with the same grouping factor. +We used the formula `y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)` to demonstrate that {{< pkg MixedModels.jl >}} amalgamates random effects terms with the same grouping factor. In practice we can use the more concise formula: `y ~ 1 + service + (1 | d) + (1 | s) + zerocorr(service | dept)`. -That is, the term `zerocorr(service | dept)` is equivalent to adding the two terms `(1 | dept) + (0 + service | dept)`, and the same as the [lme4]{.pkg} formula term `(1 + service || dept)`. +That is, the term `zerocorr(service | dept)` is equivalent to adding the two terms `(1 | dept) + (0 + service | dept)`, and the same as the {{< pkg lme4 >}} formula term `(1 + service || dept)`. -Although we only demonstrated the fitting of linear mixed models, [MixedModels.jl]{.pkg} can also fit generalized linear mixed models including logistic, probit, poisson, and negative binomial models. +Although we only demonstrated the fitting of linear mixed models, {{< pkg MixedModels.jl >}} can also fit generalized linear mixed models including logistic, probit, poisson, and negative binomial models. To fit a probit model for instance we use the following syntax: ```julia probmod = fit(MixedModel, frm, dt, Bernoulli(), ProbitLink()) ``` +\noindent where `frm` is a formula and `dt` is a dataset. -Refer to the [MixedModels.jl]{.pkg} documentation for more details. +Refer to the {{< pkg MixedModels.jl >}} documentation for more details. -## Summary and discussion {#sec-summary} +# Summary and discussion {#sec-summary} -We demonstrated how the [MixedModels.jl]{.pkg} package can be used to fit models with two distinct grouping factors (the `ml-32m` example) and three distinct grouping factors (`insteval` example). +We demonstrated how the {{< pkg MixedModels.jl >}} package can be used to fit models with two distinct grouping factors (the `ml-32m` example) and three distinct grouping factors (the `insteval` example). Fitting models with a single grouping factor or with more than three grouping factors is straightforward. -The key contribution of [MixedModels.jl]{.pkg} is the new fitting procedure based on the augmented matrix $\mbfA$ (@eq-A) and the simplified formula for the profiled loglikelihood (@eq-profiled-log-likelihood). +The key contribution of {{< pkg MixedModels.jl >}} is the new fitting procedure based on the augmented matrix $\mbfA$ (@eq-A) and the simplified formula for the profiled loglikelihood (@eq-profiled-log-likelihood). As a result of this simplification, we do not need to compute the conditional mode of the random effects ($\tilde{\mbfu}$ or $\tilde{\mbfb}$) for each evaluation of the objective function. This means we avoid roughly $O(n)$ computations at each evaluation as compared to previous implementations of mixed effects models, making the algorithm more scalable. In fact the conditional modes are computed only once at the converged values of the underlying parameter $\mbftheta$. Moreover, at each iteration we only need to use the augmented Gram matrix $\mbfA$ (@eq-A), which is a blocked, symmetric, sparse matrix of size $q + p + 1$ and the covariance parameter vector, $\mbftheta$, to form the blocked Cholesky factor, $\mbfL$. -This update is performed in-place. +This update is performed in-place in the storage for $\mbfL$. Secondly, by rearranging the random effects terms so that the grouping factor with the greatest number of random effects corresponds to the $[1,1]$ block of the Cholesky factor, we reduce the "fill-in" of the Cholesky factor. These two changes together contribute to a substantial improvement in computational efficiency. -Finally, the multiple dispatch feature of the [Julia]{.proglang} language allows us to use specialized algorithms for different kinds of matrices instead of relying on generic algorithms. +Finally, the multiple dispatch feature of the {{< proglang Julia >}} language allows us to use specialized algorithms for different kinds of matrices instead of relying on generic algorithms. This is relevant because we have both sparse and dense operations in different blocks of the Cholesky factor. By dispatching on the type of the block, we can create performant algorithms for each combination of matrix types that can occur. -## Acknowledgments {.unnumbered} +# Acknowledgments {.unnumbered} This research was supported by the Center for Interdisciplinary Research, Bielefeld (ZiF) Cooperation Group "Statistical models for psychological and linguistic data". We would like to thank Ben Bolker for helpful comments on the manuscript. -All figures were produced using the [AlgebraOfGraphics.jl]{.pkg} package with the [CairoMakie.jl]{.pkg} backend in the [Makie.jl]{.pkg} graphics system [@DanischKrumbiegel2021]. -Julius Krumbiegel and Anshul Singhvi provided invaluable help and feedback with [Makie.jl]{.pkg}. +All figures were produced using the {{< pkg AlgebraOfGraphics.jl >}} package with the {{< pkg CairoMakie.jl >}} backend in the {{< pkg Makie.jl >}} graphics system [@DanischKrumbiegel2021]. +Julius Krumbiegel and Anshul Singhvi provided invaluable help and feedback with {{< pkg Makie.jl >}}. The plots in this work would not have been possible without them. -Finally, we would like to thank the [Julia]{.proglang} community for bugfixes and answering our questions in the community support channels. +Finally, we would like to thank the {{< proglang Julia >}} community for bugfixes and answering our questions in the community support channels. -## References {.unnumbered} +# References {.unnumbered} :::{#refs} ::: - -{{< pagebreak >}} - - -\appendix -\renewcommand{\thefigure}{A\arabic{figure}} -\renewcommand{\thetable}{A\arabic{table}} -\setcounter{figure}{0} -\setcounter{table}{0} - -## Appendix {.appendix} - -### Comparative benchmarks with lme4 and glmmTMB - -We benchmarked fitting the same model as in @sec-structure using the [lme4]{.pkg} and [glmmTMB]{.pkg} packages for [R]{.proglang}. - -``` -{{< include scripts/benchmarkoutput.txt >}} -``` - -As of writing, an attempt to fit the `ml-32m` example using [glmmTMB]{.pkg} has been running for more than a week on the machine that was used for the [MixedModels.jl]{.pkg} timings and has not yet completed. -Given otherwise similar performance, we do not expect [lme4]{.pkg} to perform any better than [glmmTMB]{.pkg}. - -### Ordering of the random effects {#sec-app-re .appendix} - -To demonstrate the effect of ordering of the random effects, we fit the `insteval` model after a different ordering for the random effects. -We first recreate a copy of the `insteval` model - -```{julia} -#| eval: false -m2 = fit( - MixedModel, - @formula(y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)), - dat, - progress=false # suppress the display of a progress bar -); -``` - -To demonstrate, we must manually reorder the random effects then recreate the augmented matrix $\mbfA$ and storage for the Cholesky factor $\mbfL$ as follows: -```{julia} -#| eval: false -m2.reterms[1:2] = m2.reterms[2:-1:1]; # swap elements 1 and 2 - -A, L = MixedModels.createAL(m2.reterms, m2.Xymat); # recreate A, L - -copyto!(m2.A, A); # set A -copyto!(m2.L, L); # set L -``` - -Now we can refit the model: -```{julia} -@be refit!($m2; progress = false) seconds = 100 # very slow -``` - -This is an order of magnitude slower than just refitting the original model `m1`: -```{julia} -@be refit!($m1, progress = false) seconds = 8 -``` - -Note that the output is exactly the same as `m1` but it takes much longer to fit. -```{julia} -print(m2) -``` - -### MovieLens Runtimes {#sec-app-ml} - -The run times for the MovieLens models are summarized in @tbl-sizespeed. -The movie cutoff is the threshold on the number of ratings per movie; the user cutoff is the threshold on the number of ratings per user. The various ns are the number of ratings, users and movies in the resulting trimmed data set; the remaining columns summarize the memory and time requirements of the associated model. - - -```{=latex} -\ifXeTeX -\footnotesize -\else -\small -\fi -``` - -```{julia} -#| echo: false -#| output: asis -#| tbl-cap: Summary of results for the `ml-32m` data modeling run speeds -#| label: tbl-sizespeed -sizespeed = CSV.read("./data/sizespeed.csv", DataFrame; downcast=true); -rename!(sizespeed, - :mc => "movie cutoff", - :uc => "user cutoff", - :nratings => "ratings", - :nusers => "users", - :nmvie => "movies", - :modelsz => "model (GiB)", - :L22sz => "L[2,2] (GiB)", - :fittime => "time (s)", - :nv => "n eval", - :evtime => "time per eval (s)") -pretty_table(sizespeed; formatters = ft_round(2, [6, 7, 9, 10]), backend = Val(:markdown), show_subheader = false) -``` - -```{=latex} -\normalsize -\clearpage -``` - -### Dimensions of quantities in linear mixed models {#sec-re-dim} - -The following table is mostly a reproduction of Table 3 from @bates.maechler.etal:2015, with slight modifications. - -| Symbol | Description | `insteval` example | -|:--------------------|:--------------------------------------------------------------|:---------------------------| -| $n$ | length of the response vector $\mcY$ | $n = 73{,}421$ | -| $p$ | no. of fixed effects in the model matrix $\mbfX$ | $p = 2$ | -| $k$ | no. of random effects terms | $k = 3$ (`s`, `d`, `dept`) | -| $p_i$ | no. of fixed effects terms in the $i^\text{th}$ random effect | $(1,1,2)$ | -| $\ell_i$ | no. of levels of the $i^\text{th}$ grouping factor | $(2972, 1128, 14)$ | -| $q_i = p_i \ell_i$ | no. of columns of the term-wise model matrix | $(2972, 1128, 28)$ | -| $q = \sum_{i}q_i$ | no. of columns in the random effects model matrix $\mbfZ$ | $q = 4{,}128$ | -| $q + p + 1$ | no. of columns (and rows) of the augmented matrix $\mbfA$ | $4131$ | - - -: Sizes of model components {#tbl-re-sizes tbl-colwidths="[18,62,20]"} - -### Computing environment {#sec-app-env} - -The computing environment for all the timings, except for the models fit to the `ml-32m` data, was - -```{julia} -versioninfo() -``` - -with an accelerated BLAS library - -```{julia} -BLAS.get_config() -``` - -The versions of the Julia packages used are - -```{julia} -using Pkg; Pkg.status() -``` diff --git a/blockedcholesky/BlockedCholeskyMM.zip b/blockedcholesky/BlockedCholeskyMM.zip new file mode 100644 index 0000000..5b394b9 Binary files /dev/null and b/blockedcholesky/BlockedCholeskyMM.zip differ diff --git a/blockedcholesky/Environment.qmd b/blockedcholesky/Environment.qmd new file mode 100644 index 0000000..49b259e --- /dev/null +++ b/blockedcholesky/Environment.qmd @@ -0,0 +1,119 @@ +--- +title: "Computing environment and similar model fits in R" +format: pdf +author: + - name: Douglas Bates + email: dmbates@gmail.com + orcid: 0000-0001-8316-9503 + affiliation: + - name: University of Wisconsin--Madison + city: University of Wisconsin--Madison + state: WI + url: https://www.wisc.edu + department: Department of Statistics + - name: Phillip M. Alday + email: phillip@phillipalday.com + orcid: 0000-0002-9984-5745 + affiliation: + - name: Beacon Biosignals + city: Beacon Biosignals + url: https://beacon.bio + - name: Ajinkya H. Kokandakar + email: ajinkya@stat.wisc.edu + orcid: 0000-0001-6628-2272 + affiliations: + - name: University of Wisconsin--Madison + city: University of Wisconsin--Madison + state: WI + url: https://www.wisc.edu + department: Department of Statistics +monofont: JuliaMono +date: last-modified +number-sections: true +execute: + cache: true +engine: julia +julia: + exeflags: + - -tauto + - --project=@. +--- + +# Computing environment + +The timings in the paper, with the exception of the fits of models to the `ml-32m` data, were performed on a laptop with the following characteristics + +```{julia} +using AppleAccelerate, InteractiveUtils, LinearAlgebra, RCall +versioninfo() +``` + +The versions of the Julia packages in use were + +```{julia} +import Pkg +Pkg.status() +``` + +and the accelerated BLAS used for this computer were + +```{julia} +BLAS.get_config() +``` + +# Timing similar model fits with lme4 + +For comparison we fit the same model as in section 3.1.1 using the `lme4::lmer` and `glmmTMB::tmbmod` functions in `R` + +```{r} +dat <- lme4::InstEval +dat$service <- as.numeric(dat$service) - 1 # convert to a 0/1 numeric vector +ctrl <- lme4::lmerControl(calc.derivs = FALSE) +form <- y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) +mlmer <- lme4::lmer(form, dat, REML = FALSE, control = ctrl) +print(summary(mlmer, correlation = FALSE)) +NULL # suppress printing of the internal structure of the summary result +``` + +```{r} +tmbmod <- glmmTMB::glmmTMB(form, dat) +print(summary(tmbmod)) +NULL +``` + +These model fits are essentially the same as that from `MixedModels.jl` but they take longer + +```{r} +#| warning: false +print( + microbenchmark::microbenchmark( + lmer = lme4::lmer(form, dat, REML = FALSE, control = ctrl), + glmmTMB = glmmTMB::glmmTMB(form, dat), + times=6L, + unit="seconds" + ) +) +NULL +``` + +Recall from section 3.1.1 of the paper that the `MixedModels.jl` fit of this model took just over one second on the same laptop computer. + +We attempted to fit the `ml-32m` example using `glmmTMB` on a desktop computer with 64 GiB of memory but shut it down after it had been running for more than a week on the machine without converging. +Given otherwise similar performance, we do not expect `lme4` to perform any better than `glmmTMB`. + +# Versions of R and the R packages that were used + +```{r} +print(R.version.string) +print(packageVersion("lme4")) +print(packageVersion("glmmTMB")) +NULL +``` + +# Environment for the ml-32m model fits + +The environment for the ml-32m model fits is shown in the `config.txt` file. The Julia source code to download the `ml-32m.zip` file, extract the `ratings` table, create `moviesnratings.arrow` and `usersnratings.arrow`, followed by `ratings.arrow` is `ml32-download.jl`. + +The Julia source file to run the `ml-32m` fits is `ml32-timefit.jl`. +Before you start these fits remember that they take up to a day on a powerful server computer. +Don't do this frivolously. diff --git a/blockedcholesky/Makefile b/blockedcholesky/Makefile index 9363d41..01b0970 100644 --- a/blockedcholesky/Makefile +++ b/blockedcholesky/Makefile @@ -7,7 +7,7 @@ preview: quarto preview $(DOCNAME).qmd render: - quarto render $(DOCNAME).qmd --to jss-pdf+jss + quarto render $(DOCNAME).qmd --to jds-pdf render-arxiv: sed -i'.bak' -e '/city:/d' $(DOCNAME).qmd diff --git a/blockedcholesky/Project.toml b/blockedcholesky/Project.toml index b9b39ba..51cd1ac 100644 --- a/blockedcholesky/Project.toml +++ b/blockedcholesky/Project.toml @@ -6,21 +6,26 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7" MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" +PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +RCall = "6f49c342-dc21-5d91-9882-a32aef131414" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [compat] -AlgebraOfGraphics = "0.10.4" +AlgebraOfGraphics = "0.11,0.12" Arrow = "2.8" CSV = "0.10" -CairoMakie = "0.13.4" +CairoMakie = "0.13,0.14,0.15" Chairmarks = "1.3" DataFrames = "1.7" -MixedModels = "4.34" -PrettyTables = "2.4" +MixedModels = "4,5" +PrettyTables = "2,3" Printf = "1" -julia = "1.10" +julia = "1.12" diff --git a/blockedcholesky/README.md b/blockedcholesky/README.md index dabaae6..58059e1 100644 --- a/blockedcholesky/README.md +++ b/blockedcholesky/README.md @@ -3,26 +3,31 @@ ## Steps -- Install (or update) [`quarto`](https://quarto.org/) version 1.7 or greater (v1.6 may work as well, basically use any version of `quarto` that supports `engine: julia`) +- Install the current Julia release version according to the instructions at https://julialang.org/downloads +- Install (or update) [`quarto`](https://quarto.org/) version 1.9.0 or greater - Install `jupyter-cache` using the instructions from the quarto Julia docs [here](https://quarto.org/docs/computations/julia.html#jupyter-cache). - Install the `JuliaMono` font, see instructions at: https://juliamono.netlify.app -- Install required `Julia` dependencies by running: - ```bash - $ julia --project=@. -e "import Pkg; Pkg.instantiate()" - ``` -- Generate the JSS formatted preview using: +- Install the [`JDS extensions`](https://github.com/wenjie2wang/jds.qmd) for [`quarto`](https://quarto.org) +- Install the required `Julia` dependencies by running: ```bash - $ make render + julia --project=@. -e "import Pkg; Pkg.instantiate()" ``` -- Generate the arXiv formatted version using: +- Generate the JDL formatted preview using: ```bash - $ make render-arxiv + make render ``` - **Note:** this command re-runs the `R` benchmark script which can take a long time. - **Note:** For the arxiv version, you need to have `cbfonts-fd` and `cbfonts` packages installed in your latex distribution - If you get errors due to caching, try using the `--no-cache` option for quarto: ```bash - quarto render BlockedCholeskyMM.qmd --no-cache # for JSS version render - quarto render BlockedCholeskyMM.qmd --no-cache --to arxiv-pdf+arxiv # for arxiv version - + quarto render BlockedCholeskyMM.qmd --no-cache --to jds-pdf # for JDS version render + ``` +- To create a Julia script from the Quarto file (`BlockedCholeskyMM.qmd`) convert the .qmd file to a Jupyter notebook and apply `jupyter nbconvert` to create a script. + ```bash + quarto convert BlockedCholeskyMM.qmd + jupyter nbconvert --no-prompt --to script BlockedCholeskyMM.ipynb + mv BlockedCholeskyMM.txt BlockedCholeskyMM.jl ``` +- Alternatively, `Quarto mode` for editors such as [VSCode](https://code.visualstudio.com) or [Positron](https://positron.posit.co) provide the ability to interactively evaluate code chunks from a Quarto document. + +- Comparison fits in R using the `lme4` and `glmmTMB` packages are shown in `Environment.pdf` generated from `Environment.qmd` + +- Because the ml-32m dataset cannot be redistributed, reproduction of the results in Table 4, require a separate script, `scripts/ml32-download.jl` to download the data, create the ratings table, and save it in the [Arrow](https://arrow.apache.org) format. The evaluation of the models is performed by `ml32-timefit.jl` **but** be aware that fitting very large models like this requires a computer with a large amount of memory (64 GiB is recommended) and can take a long time. It took roughly 18 hours for all the model fits on a cloud compute instance with the characteristics shown in `config.txt`. diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/_extension.yml b/blockedcholesky/_extensions/wenjie2wang/jds/_extension.yml new file mode 100644 index 0000000..0b85695 --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/_extension.yml @@ -0,0 +1,27 @@ +title: Journal of Data Science +author: Wenjie Wang and Jun Yan +version: 0.1.1 +quarto-required: ">= 1.4.0" +contributes: + formats: + common: + shortcodes: + - shortcodes.lua + pdf: + documentclass: jdsart + classoption: + - letterpaper + - inpress + template: template.tex + cite-method: natbib + biblio-style: jds + biblio-config: false + fig-width: 6.5 + fig-height: 4.5 + format-resources: + - jdsart.cls + - jds.bst + html: + number-sections: true + toc: true + csl: apa.csl diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/apa.csl b/blockedcholesky/_extensions/wenjie2wang/jds/apa.csl new file mode 100644 index 0000000..84632ab --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/apa.csl @@ -0,0 +1,2273 @@ + + diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/jds.bst b/blockedcholesky/_extensions/wenjie2wang/jds/jds.bst new file mode 100644 index 0000000..92049fd --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/jds.bst @@ -0,0 +1,1537 @@ +%% +%% This is file `jds.bst', +%% generated with the docstrip utility. +%% +%% The original source files were: +%% +%% merlin.mbs (with options: `ay,nat,nm-rvx,jnrlst,nmlm,x6,m6,keyxyr,dt-beg,yr-par,note-yr,thtit-a,trnum-it,volp-sp,num-xser,edparxc,bkedcap,edby-par,edbyy,in-col,ppx,ed,xedn,xand,nfss,') +%% ---------------------------------------- +%% *** Bibstyle for the Journal of Data Science *** +%% +%% Copyright (C) 2020 Journal of Data Science Editorial Office +%% +%% This work may be distributed and/or modified under the +%% conditions of the LaTeX Project Public License, either version 1.3 +%% of this license or (at your option) any later version. +%% The latest version of this license is in +%% +%% and version 1.3 or later is part of all distributions of LaTeX +%% version 2005/12/01 or later. +%% This work has the LPPL maintenance status `maintained'. +%% +%% The Current Maintainer of this work is Jun Yan. +%% This work consists of the file jds.dbj and the derived file jds.bst. +%% +%% Copyright 1994-2011 Patrick W Daly + % =============================================================== + % IMPORTANT NOTICE: + % This bibliographic style (bst) file has been generated from one or + % more master bibliographic style (mbs) files, listed above. + % + % This generated file can be redistributed and/or modified under the terms + % of the LaTeX Project Public License Distributed from CTAN + % archives in directory macros/latex/base/lppl.txt; either + % version 1 of the License, or any later version. + % =============================================================== + % Name and version information of the main mbs file: + % \ProvidesFile{merlin.mbs}[2011/11/18 4.33 (PWD, AO, DPC)] + % For use with BibTeX version 0.99a or later + %------------------------------------------------------------------- + % This bibliography style file is intended for texts in ENGLISH + % This is an author-year citation style bibliography. As such, it is + % non-standard LaTeX, and requires a special package file to function properly. + % Such a package is natbib.sty by Patrick W. Daly + % The form of the \bibitem entries is + % \bibitem[Jones et al.(1990)]{key}... + % \bibitem[Jones et al.(1990)Jones, Baker, and Smith]{key}... + % The essential feature is that the label (the part in brackets) consists + % of the author names, as they should appear in the citation, with the year + % in parentheses following. There must be no space before the opening + % parenthesis! + % With natbib v5.3, a full list of authors may also follow the year. + % In natbib.sty, it is possible to define the type of enclosures that is + % really wanted (brackets or parentheses), but in either case, there must + % be parentheses in the label. + % The \cite command functions as follows: + % \citet{key} ==>> Jones et al. (1990) + % \citet*{key} ==>> Jones, Baker, and Smith (1990) + % \citep{key} ==>> (Jones et al., 1990) + % \citep*{key} ==>> (Jones, Baker, and Smith, 1990) + % \citep[chap. 2]{key} ==>> (Jones et al., 1990, chap. 2) + % \citep[e.g.][]{key} ==>> (e.g. Jones et al., 1990) + % \citep[e.g.][p. 32]{key} ==>> (e.g. Jones et al., 1990, p. 32) + % \citeauthor{key} ==>> Jones et al. + % \citeauthor*{key} ==>> Jones, Baker, and Smith + % \citeyear{key} ==>> 1990 + %--------------------------------------------------------------------- + +ENTRY + { address + author + booktitle + chapter + edition + editor + eid + howpublished + institution + journal + key + month + note + number + organization + pages + publisher + school + series + title + type + volume + year + } + {} + { label extra.label sort.label short.list } +INTEGERS { output.state before.all mid.sentence after.sentence after.block } +FUNCTION {init.state.consts} +{ #0 'before.all := + #1 'mid.sentence := + #2 'after.sentence := + #3 'after.block := +} +STRINGS { s t} +FUNCTION {output.nonnull} +{ 's := + output.state mid.sentence = + { ", " * write$ } + { output.state after.block = + { add.period$ write$ + newline$ + "\newblock " write$ + } + { output.state before.all = + 'write$ + { add.period$ " " * write$ } + if$ + } + if$ + mid.sentence 'output.state := + } + if$ + s +} +FUNCTION {output} +{ duplicate$ empty$ + 'pop$ + 'output.nonnull + if$ +} +FUNCTION {output.check} +{ 't := + duplicate$ empty$ + { pop$ "empty " t * " in " * cite$ * warning$ } + 'output.nonnull + if$ +} +FUNCTION {fin.entry} +{ add.period$ + write$ + newline$ +} + +FUNCTION {new.block} +{ output.state before.all = + 'skip$ + { after.block 'output.state := } + if$ +} +FUNCTION {new.sentence} +{ output.state after.block = + 'skip$ + { output.state before.all = + 'skip$ + { after.sentence 'output.state := } + if$ + } + if$ +} +FUNCTION {add.blank} +{ " " * before.all 'output.state := +} + +FUNCTION {date.block} +{ + new.block +} + +FUNCTION {not} +{ { #0 } + { #1 } + if$ +} +FUNCTION {and} +{ 'skip$ + { pop$ #0 } + if$ +} +FUNCTION {or} +{ { pop$ #1 } + 'skip$ + if$ +} +STRINGS {z} + +FUNCTION {remove.dots} +{ 'z := + "" + { z empty$ not } + { z #1 #2 substring$ + duplicate$ "\." = + { z #3 global.max$ substring$ 'z := * } + { pop$ + z #1 #1 substring$ + z #2 global.max$ substring$ 'z := + duplicate$ "." = 'pop$ + { * } + if$ + } + if$ + } + while$ +} +FUNCTION {new.block.checkb} +{ empty$ + swap$ empty$ + and + 'skip$ + 'new.block + if$ +} +FUNCTION {field.or.null} +{ duplicate$ empty$ + { pop$ "" } + 'skip$ + if$ +} +FUNCTION {emphasize} +{ duplicate$ empty$ + { pop$ "" } + { "\emph{" swap$ * "}" * } + if$ +} +FUNCTION {tie.or.space.prefix} +{ duplicate$ text.length$ #3 < + { "~" } + { " " } + if$ + swap$ +} + +FUNCTION {capitalize} +{ "u" change.case$ "t" change.case$ } + +FUNCTION {space.word} +{ " " swap$ * " " * } + % Here are the language-specific definitions for explicit words. + % Each function has a name bbl.xxx where xxx is the English word. + % The language selected here is ENGLISH +FUNCTION {bbl.and} +{ "and"} + +FUNCTION {bbl.etal} +{ "et~al." } + +FUNCTION {bbl.editors} +{ "eds." } + +FUNCTION {bbl.editor} +{ "ed." } + +FUNCTION {bbl.edby} +{ "edited by" } + +FUNCTION {bbl.edition} +{ "edition" } + +FUNCTION {bbl.volume} +{ "volume" } + +FUNCTION {bbl.of} +{ "of" } + +FUNCTION {bbl.number} +{ "number" } + +FUNCTION {bbl.nr} +{ "no." } + +FUNCTION {bbl.in} +{ "in" } + +FUNCTION {bbl.pages} +{ "" } + +FUNCTION {bbl.page} +{ "" } + +FUNCTION {bbl.chapter} +{ "chapter" } + +FUNCTION {bbl.techrep} +{ "Technical Report" } + +FUNCTION {bbl.mthesis} +{ "Master's thesis" } + +FUNCTION {bbl.phdthesis} +{ "Ph.D. thesis" } + +MACRO {jan} {"January"} + +MACRO {feb} {"February"} + +MACRO {mar} {"March"} + +MACRO {apr} {"April"} + +MACRO {may} {"May"} + +MACRO {jun} {"June"} + +MACRO {jul} {"July"} + +MACRO {aug} {"August"} + +MACRO {sep} {"September"} + +MACRO {oct} {"October"} + +MACRO {nov} {"November"} + +MACRO {dec} {"December"} + +MACRO {acmcs} {"ACM Computing Surveys"} + +MACRO {acta} {"Acta Informatica"} + +MACRO {cacm} {"Communications of the ACM"} + +MACRO {ibmjrd} {"IBM Journal of Research and Development"} + +MACRO {ibmsj} {"IBM Systems Journal"} + +MACRO {ieeese} {"IEEE Transactions on Software Engineering"} + +MACRO {ieeetc} {"IEEE Transactions on Computers"} + +MACRO {ieeetcad} + {"IEEE Transactions on Computer-Aided Design of Integrated Circuits"} + +MACRO {ipl} {"Information Processing Letters"} + +MACRO {jacm} {"Journal of the ACM"} + +MACRO {jcss} {"Journal of Computer and System Sciences"} + +MACRO {scp} {"Science of Computer Programming"} + +MACRO {sicomp} {"SIAM Journal on Computing"} + +MACRO {tocs} {"ACM Transactions on Computer Systems"} + +MACRO {tods} {"ACM Transactions on Database Systems"} + +MACRO {tog} {"ACM Transactions on Graphics"} + +MACRO {toms} {"ACM Transactions on Mathematical Software"} + +MACRO {toois} {"ACM Transactions on Office Information Systems"} + +MACRO {toplas} {"ACM Transactions on Programming Languages and Systems"} + +MACRO {tcs} {"Theoretical Computer Science"} +FUNCTION {bibinfo.check} +{ swap$ + duplicate$ missing$ + { + pop$ pop$ + "" + } + { duplicate$ empty$ + { + swap$ pop$ + } + { swap$ + pop$ + } + if$ + } + if$ +} +FUNCTION {bibinfo.warn} +{ swap$ + duplicate$ missing$ + { + swap$ "missing " swap$ * " in " * cite$ * warning$ pop$ + "" + } + { duplicate$ empty$ + { + swap$ "empty " swap$ * " in " * cite$ * warning$ + } + { swap$ + pop$ + } + if$ + } + if$ +} +INTEGERS { nameptr namesleft numnames } + + +STRINGS { bibinfo} + +FUNCTION {format.names} +{ 'bibinfo := + duplicate$ empty$ 'skip$ { + 's := + "" 't := + #1 'nameptr := + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr + "{vv~}{ll}{ f{}}{ jj}" + format.name$ + remove.dots + bibinfo bibinfo.check + 't := + nameptr #1 > + { + nameptr #6 + #1 + = + numnames #6 + > and + { "others" 't := + #1 'namesleft := } + 'skip$ + if$ + namesleft #1 > + { ", " * t * } + { + s nameptr "{ll}" format.name$ duplicate$ "others" = + { 't := } + { pop$ } + if$ + "," * + t "others" = + { + " " * bbl.etal * + } + { " " * t * } + if$ + } + if$ + } + 't + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ + } if$ +} +FUNCTION {format.names.ed} +{ + 'bibinfo := + duplicate$ empty$ 'skip$ { + 's := + "" 't := + #1 'nameptr := + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr + "{f{}~}{vv~}{ll}{ jj}" + format.name$ + remove.dots + bibinfo bibinfo.check + 't := + nameptr #1 > + { + namesleft #1 > + { ", " * t * } + { + s nameptr "{ll}" format.name$ duplicate$ "others" = + { 't := } + { pop$ } + if$ + "," * + t "others" = + { + + " " * bbl.etal * + } + { " " * t * } + if$ + } + if$ + } + 't + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ + } if$ +} +FUNCTION {format.key} +{ empty$ + { key field.or.null } + { "" } + if$ +} + +FUNCTION {format.authors} +{ author "author" format.names +} +FUNCTION {get.bbl.editor} +{ editor num.names$ #1 > 'bbl.editors 'bbl.editor if$ } + +FUNCTION {format.editors} +{ editor "editor" format.names duplicate$ empty$ 'skip$ + { + " " * + get.bbl.editor + capitalize + "(" swap$ * ")" * + * + } + if$ +} +FUNCTION {format.note} +{ + note empty$ + { "" } + { note #1 #1 substring$ + duplicate$ "{" = + 'skip$ + { output.state mid.sentence = + { "l" } + { "u" } + if$ + change.case$ + } + if$ + note #2 global.max$ substring$ * "note" bibinfo.check + } + if$ +} + +FUNCTION {format.title} +{ title + duplicate$ empty$ 'skip$ + { "t" change.case$ } + if$ + "title" bibinfo.check +} +FUNCTION {format.full.names} +{'s := + "" 't := + #1 'nameptr := + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr + "{vv~}{ll}" format.name$ + 't := + nameptr #1 > + { + nameptr #6 + #1 + = + numnames #6 + > and + { "others" 't := + #1 'namesleft := } + 'skip$ + if$ + namesleft #1 > + { ", " * t * } + { + s nameptr "{ll}" format.name$ duplicate$ "others" = + { 't := } + { pop$ } + if$ + t "others" = + { + " " * bbl.etal * + } + { + numnames #2 > + { "," * } + 'skip$ + if$ + bbl.and + space.word * t * + } + if$ + } + if$ + } + 't + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ +} + +FUNCTION {author.editor.key.full} +{ author empty$ + { editor empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { editor format.full.names } + if$ + } + { author format.full.names } + if$ +} + +FUNCTION {author.key.full} +{ author empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { author format.full.names } + if$ +} + +FUNCTION {editor.key.full} +{ editor empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { editor format.full.names } + if$ +} + +FUNCTION {make.full.names} +{ type$ "book" = + type$ "inbook" = + or + 'author.editor.key.full + { type$ "proceedings" = + 'editor.key.full + 'author.key.full + if$ + } + if$ +} + +FUNCTION {output.bibitem} +{ newline$ + "\bibitem[{" write$ + label write$ + ")" make.full.names duplicate$ short.list = + { pop$ } + { * } + if$ + "}]{" * write$ + cite$ write$ + "}" write$ + newline$ + "" + before.all 'output.state := +} + +FUNCTION {n.dashify} +{ + 't := + "" + { t empty$ not } + { t #1 #1 substring$ "-" = + { t #1 #2 substring$ "--" = not + { "--" * + t #2 global.max$ substring$ 't := + } + { { t #1 #1 substring$ "-" = } + { "-" * + t #2 global.max$ substring$ 't := + } + while$ + } + if$ + } + { t #1 #1 substring$ * + t #2 global.max$ substring$ 't := + } + if$ + } + while$ +} + +FUNCTION {word.in} +{ bbl.in capitalize + ":" * + " " * } + +FUNCTION {format.date} +{ year "year" bibinfo.check duplicate$ empty$ + { + "empty year in " cite$ * "; set to ????" * warning$ + pop$ "????" + } + 'skip$ + if$ + extra.label * + before.all 'output.state := + " (" swap$ * ")" * +} +FUNCTION {format.btitle} +{ title "title" bibinfo.check + duplicate$ empty$ 'skip$ + { + emphasize + } + if$ +} +FUNCTION {either.or.check} +{ empty$ + 'pop$ + { "can't use both " swap$ * " fields in " * cite$ * warning$ } + if$ +} +FUNCTION {format.bvolume} +{ volume empty$ + { "" } + { bbl.volume volume tie.or.space.prefix + "volume" bibinfo.check * * + series "series" bibinfo.check + duplicate$ empty$ 'pop$ + { swap$ bbl.of space.word * swap$ + emphasize * } + if$ + "volume and number" number either.or.check + } + if$ +} +FUNCTION {format.number.series} +{ volume empty$ + { number empty$ + { series field.or.null } + { series empty$ + { number "number" bibinfo.check } + { output.state mid.sentence = + { bbl.number } + { bbl.number capitalize } + if$ + number tie.or.space.prefix "number" bibinfo.check * * + bbl.in space.word * + series "series" bibinfo.check * + } + if$ + } + if$ + } + { "" } + if$ +} + +FUNCTION {format.edition} +{ edition duplicate$ empty$ 'skip$ + { + output.state mid.sentence = + { "l" } + { "t" } + if$ change.case$ + "edition" bibinfo.check + " " * bbl.edition * + } + if$ +} +INTEGERS { multiresult } +FUNCTION {multi.page.check} +{ 't := + #0 'multiresult := + { multiresult not + t empty$ not + and + } + { t #1 #1 substring$ + duplicate$ "-" = + swap$ duplicate$ "," = + swap$ "+" = + or or + { #1 'multiresult := } + { t #2 global.max$ substring$ 't := } + if$ + } + while$ + multiresult +} +FUNCTION {format.pages} +{ pages duplicate$ empty$ 'skip$ + { duplicate$ multi.page.check + { + n.dashify + } + { + } + if$ + "pages" bibinfo.check + } + if$ +} +FUNCTION {format.journal.pages} +{ pages duplicate$ empty$ 'pop$ + { swap$ duplicate$ empty$ + { pop$ pop$ format.pages } + { + ": " * + swap$ + n.dashify + "pages" bibinfo.check + * + } + if$ + } + if$ +} +FUNCTION {format.journal.eid} +{ eid "eid" bibinfo.check + duplicate$ empty$ 'pop$ + { swap$ duplicate$ empty$ 'skip$ + { + ": " * + } + if$ + swap$ * + } + if$ +} +FUNCTION {format.vol.num.pages} +{ volume field.or.null + duplicate$ empty$ 'skip$ + { + "volume" bibinfo.check + } + if$ + number "number" bibinfo.check duplicate$ empty$ 'skip$ + { + swap$ duplicate$ empty$ + { "there's a number but no volume in " cite$ * warning$ } + 'skip$ + if$ + swap$ + "(" swap$ * ")" * + } + if$ * + eid empty$ + { format.journal.pages } + { format.journal.eid } + if$ +} + +FUNCTION {format.chapter.pages} +{ chapter empty$ + 'format.pages + { type empty$ + { bbl.chapter } + { type "l" change.case$ + "type" bibinfo.check + } + if$ + chapter tie.or.space.prefix + "chapter" bibinfo.check + * * + pages empty$ + 'skip$ + { ", " * format.pages * } + if$ + } + if$ +} + +FUNCTION {format.booktitle} +{ + booktitle "booktitle" bibinfo.check + emphasize +} +FUNCTION {format.in.ed.booktitle} +{ format.booktitle duplicate$ empty$ 'skip$ + { + editor "editor" format.names.ed duplicate$ empty$ 'pop$ + { + get.bbl.editor + swap$ "," * + " " * swap$ * + "(" swap$ * ")" * + swap$ + " " * swap$ + * } + if$ + word.in swap$ * + } + if$ +} +FUNCTION {format.thesis.type} +{ type duplicate$ empty$ + 'pop$ + { swap$ pop$ + "t" change.case$ "type" bibinfo.check + } + if$ +} +FUNCTION {format.tr.number} +{ number "number" bibinfo.check + type duplicate$ empty$ + { pop$ bbl.techrep } + 'skip$ + if$ + "type" bibinfo.check + swap$ duplicate$ empty$ + { pop$ "t" change.case$ } + { tie.or.space.prefix * * } + if$ +} +FUNCTION {format.article.crossref} +{ + word.in + " \cite{" * crossref * "}" * +} +FUNCTION {format.book.crossref} +{ volume duplicate$ empty$ + { "empty volume in " cite$ * "'s crossref of " * crossref * warning$ + pop$ word.in + } + { bbl.volume + capitalize + swap$ tie.or.space.prefix "volume" bibinfo.check * * bbl.of space.word * + } + if$ + " \cite{" * crossref * "}" * +} +FUNCTION {format.incoll.inproc.crossref} +{ + word.in + " \cite{" * crossref * "}" * +} +FUNCTION {format.org.or.pub} +{ 't := + "" + address empty$ t empty$ and + 'skip$ + { + t empty$ + { address "address" bibinfo.check * + } + { t * + address empty$ + 'skip$ + { ", " * address "address" bibinfo.check * } + if$ + } + if$ + } + if$ +} +FUNCTION {format.publisher.address} +{ publisher "publisher" bibinfo.warn format.org.or.pub +} + +FUNCTION {format.organization.address} +{ organization "organization" bibinfo.check format.org.or.pub +} + +FUNCTION {article} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title "title" output.check + new.block + crossref missing$ + { + journal + "journal" bibinfo.check + emphasize + "journal" output.check + format.vol.num.pages output + } + { format.article.crossref output.nonnull + format.pages output + } + if$ + new.block + format.note output + fin.entry +} +FUNCTION {book} +{ output.bibitem + author empty$ + { format.editors "author and editor" output.check + editor format.key output + add.blank + } + { format.authors output.nonnull + crossref missing$ + { "author and editor" editor either.or.check } + 'skip$ + if$ + } + if$ + format.date "year" output.check + date.block + format.btitle "title" output.check + crossref missing$ + { format.bvolume output + new.block + format.number.series output + new.sentence + format.publisher.address output + } + { + new.block + format.book.crossref output.nonnull + } + if$ + format.edition output + new.block + format.note output + fin.entry +} +FUNCTION {booklet} +{ output.bibitem + format.authors output + author format.key output + format.date "year" output.check + date.block + format.title "title" output.check + new.block + howpublished "howpublished" bibinfo.check output + address "address" bibinfo.check output + new.block + format.note output + fin.entry +} + +FUNCTION {inbook} +{ output.bibitem + author empty$ + { format.editors "author and editor" output.check + editor format.key output + } + { format.authors output.nonnull + crossref missing$ + { "author and editor" editor either.or.check } + 'skip$ + if$ + } + if$ + format.date "year" output.check + date.block + format.btitle "title" output.check + crossref missing$ + { + format.bvolume output + format.chapter.pages "chapter and pages" output.check + new.block + format.number.series output + new.sentence + format.publisher.address output + } + { + format.chapter.pages "chapter and pages" output.check + new.block + format.book.crossref output.nonnull + } + if$ + format.edition output + new.block + format.note output + fin.entry +} + +FUNCTION {incollection} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title "title" output.check + new.block + crossref missing$ + { format.in.ed.booktitle "booktitle" output.check + format.bvolume output + format.number.series output + format.chapter.pages output + new.sentence + format.publisher.address output + format.edition output + } + { format.incoll.inproc.crossref output.nonnull + format.chapter.pages output + } + if$ + new.block + format.note output + fin.entry +} +FUNCTION {inproceedings} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title "title" output.check + new.block + crossref missing$ + { format.in.ed.booktitle "booktitle" output.check + format.bvolume output + format.number.series output + format.pages output + new.sentence + publisher empty$ + { format.organization.address output } + { organization "organization" bibinfo.check output + format.publisher.address output + } + if$ + } + { format.incoll.inproc.crossref output.nonnull + format.pages output + } + if$ + new.block + format.note output + fin.entry +} +FUNCTION {conference} { inproceedings } +FUNCTION {manual} +{ output.bibitem + format.authors output + author format.key output + format.date "year" output.check + date.block + format.btitle "title" output.check + organization address new.block.checkb + organization "organization" bibinfo.check output + address "address" bibinfo.check output + format.edition output + new.block + format.note output + fin.entry +} + +FUNCTION {mastersthesis} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title + "title" output.check + new.block + bbl.mthesis format.thesis.type output.nonnull + school "school" bibinfo.warn output + address "address" bibinfo.check output + new.block + format.note output + fin.entry +} + +FUNCTION {misc} +{ output.bibitem + format.authors output + author format.key output + format.date "year" output.check + date.block + format.title output + new.block + howpublished "howpublished" bibinfo.check output + new.block + format.note output + fin.entry +} +FUNCTION {phdthesis} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title + "title" output.check + new.block + bbl.phdthesis format.thesis.type output.nonnull + school "school" bibinfo.warn output + address "address" bibinfo.check output + new.block + format.note output + fin.entry +} + +FUNCTION {proceedings} +{ output.bibitem + format.editors output + editor format.key output + format.date "year" output.check + date.block + format.btitle "title" output.check + format.bvolume output + format.number.series output + new.sentence + publisher empty$ + { format.organization.address output } + { organization "organization" bibinfo.check output + format.publisher.address output + } + if$ + new.block + format.note output + fin.entry +} + +FUNCTION {techreport} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title + "title" output.check + new.block + format.tr.number emphasize output.nonnull + institution "institution" bibinfo.warn output + address "address" bibinfo.check output + new.block + format.note output + fin.entry +} + +FUNCTION {unpublished} +{ output.bibitem + format.authors "author" output.check + author format.key output + format.date "year" output.check + date.block + format.title "title" output.check + new.block + format.note "note" output.check + fin.entry +} + +FUNCTION {default.type} { misc } +READ +FUNCTION {sortify} +{ purify$ + "l" change.case$ +} +INTEGERS { len } +FUNCTION {chop.word} +{ 's := + 'len := + s #1 len substring$ = + { s len #1 + global.max$ substring$ } + 's + if$ +} +FUNCTION {format.lab.names} +{ 's := + "" 't := + s #1 "{vv~}{ll}" format.name$ + s num.names$ duplicate$ + #2 > + { pop$ + " " * bbl.etal * + } + { #2 < + 'skip$ + { s #2 "{ff }{vv }{ll}{ jj}" format.name$ "others" = + { + " " * bbl.etal * + } + { bbl.and space.word * s #2 "{vv~}{ll}" format.name$ + * } + if$ + } + if$ + } + if$ +} + +FUNCTION {author.key.label} +{ author empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { author format.lab.names } + if$ +} + +FUNCTION {author.editor.key.label} +{ author empty$ + { editor empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { editor format.lab.names } + if$ + } + { author format.lab.names } + if$ +} + +FUNCTION {editor.key.label} +{ editor empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { editor format.lab.names } + if$ +} + +FUNCTION {calc.short.authors} +{ type$ "book" = + type$ "inbook" = + or + 'author.editor.key.label + { type$ "proceedings" = + 'editor.key.label + 'author.key.label + if$ + } + if$ + 'short.list := +} + +FUNCTION {calc.label} +{ calc.short.authors + short.list + "(" + * + year duplicate$ empty$ + short.list key field.or.null = or + { pop$ "" } + 'skip$ + if$ + * + 'label := +} + +FUNCTION {sort.format.names} +{ 's := + #1 'nameptr := + "" + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr + "{vv{ } }{ll{ }}{ f{ }}{ jj{ }}" + format.name$ 't := + nameptr #1 > + { + nameptr #6 + #1 + = + numnames #6 + > and + { "others" 't := + #1 'namesleft := } + 'skip$ + if$ + " " * + namesleft #1 = t "others" = and + { "zzzzz" 't := } + 'skip$ + if$ + t sortify * + } + { t sortify * } + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ +} + +FUNCTION {sort.format.title} +{ 't := + "A " #2 + "An " #3 + "The " #4 t chop.word + chop.word + chop.word + sortify + #1 global.max$ substring$ +} +FUNCTION {author.sort} +{ author empty$ + { key empty$ + { "to sort, need author or key in " cite$ * warning$ + "" + } + { key sortify } + if$ + } + { author sort.format.names } + if$ +} +FUNCTION {author.editor.sort} +{ author empty$ + { editor empty$ + { key empty$ + { "to sort, need author, editor, or key in " cite$ * warning$ + "" + } + { key sortify } + if$ + } + { editor sort.format.names } + if$ + } + { author sort.format.names } + if$ +} +FUNCTION {editor.sort} +{ editor empty$ + { key empty$ + { "to sort, need editor or key in " cite$ * warning$ + "" + } + { key sortify } + if$ + } + { editor sort.format.names } + if$ +} +FUNCTION {presort} +{ calc.label + label sortify + " " + * + type$ "book" = + type$ "inbook" = + or + 'author.editor.sort + { type$ "proceedings" = + 'editor.sort + 'author.sort + if$ + } + if$ + #1 entry.max$ substring$ + 'sort.label := + sort.label + * + " " + * + title field.or.null + sort.format.title + * + #1 entry.max$ substring$ + 'sort.key$ := +} + +ITERATE {presort} +SORT +STRINGS { last.label next.extra } +INTEGERS { last.extra.num last.extra.num.extended last.extra.num.blank number.label } +FUNCTION {initialize.extra.label.stuff} +{ #0 int.to.chr$ 'last.label := + "" 'next.extra := + #0 'last.extra.num := + "a" chr.to.int$ #1 - 'last.extra.num.blank := + last.extra.num.blank 'last.extra.num.extended := + #0 'number.label := +} +FUNCTION {forward.pass} +{ last.label label = + { last.extra.num #1 + 'last.extra.num := + last.extra.num "z" chr.to.int$ > + { "a" chr.to.int$ 'last.extra.num := + last.extra.num.extended #1 + 'last.extra.num.extended := + } + 'skip$ + if$ + last.extra.num.extended last.extra.num.blank > + { last.extra.num.extended int.to.chr$ + last.extra.num int.to.chr$ + * 'extra.label := } + { last.extra.num int.to.chr$ 'extra.label := } + if$ + } + { "a" chr.to.int$ 'last.extra.num := + "" 'extra.label := + label 'last.label := + } + if$ + number.label #1 + 'number.label := +} +FUNCTION {reverse.pass} +{ next.extra "b" = + { "a" 'extra.label := } + 'skip$ + if$ + extra.label 'next.extra := + extra.label + duplicate$ empty$ + 'skip$ + { "{\natexlab{" swap$ * "}}" * } + if$ + 'extra.label := + label extra.label * 'label := +} +EXECUTE {initialize.extra.label.stuff} +ITERATE {forward.pass} +REVERSE {reverse.pass} +FUNCTION {bib.sort.order} +{ sort.label + " " + * + year field.or.null sortify + * + " " + * + title field.or.null + sort.format.title + * + #1 entry.max$ substring$ + 'sort.key$ := +} +ITERATE {bib.sort.order} +SORT +FUNCTION {begin.bib} +{ preamble$ empty$ + 'skip$ + { preamble$ write$ newline$ } + if$ + "\begin{thebibliography}{" number.label int.to.str$ * "}" * + write$ newline$ + "\providecommand{\natexlab}[1]{#1}" + write$ newline$ +} +EXECUTE {begin.bib} +EXECUTE {init.state.consts} +ITERATE {call.type$} +FUNCTION {end.bib} +{ newline$ + "\end{thebibliography}" write$ newline$ +} +EXECUTE {end.bib} +%% End of customized bst file +%% +%% End of file `jds.bst'. diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/jdsart.cls b/blockedcholesky/_extensions/wenjie2wang/jds/jdsart.cls new file mode 100644 index 0000000..2976658 --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/jdsart.cls @@ -0,0 +1,1903 @@ +%% This is file `jdsart.cls'. +%% +%% LaTeX 2e class file for articles in jds journal published by VTeX. +%% +%% Macros written by VTeX Ltd., Lithuania +%% Please submit bugs or your comments to latex-support@vtex.lt +%% +%% Copyright 2020 by VTEX UAB +%% +%% You are free to use this class file as you see fit, provided +%% that you do not make changes to the file. +%% If you DO make changes, you are required to rename this file. +%% +%% It may be distributed under the terms of the LaTeX Project Public +%% License, as described in lppl.txt in the base LaTeX distribution. +%% Either version 1.3 or (at your option) any later version. +%% +\NeedsTeXFormat{LaTeX2e} +\ProvidesClass{jdsart}[2025/09/03 v0.25 Journal of Data Science class] + +\def\fmt@name{jdsart} + +\RequirePackage{etoolbox} + +\newif\if@review +\newif\if@inpress \@inpresstrue +\newif\if@discussion +\newif\if@chinese +\newif\if@linenumber + +% only one of review and inpress can be true +\DeclareOption{review}{\@reviewtrue \@inpressfalse \@linenumbertrue} +\DeclareOption{inpress}{\@reviewfalse \@inpresstrue} +\DeclareOption{discussion}{\@discussiontrue} +\DeclareOption{chinese}{\@chinesetrue \@discussiontrue} +\DeclareOption{linenumber}{\@linenumbertrue} +\DeclareOption{nolinenumber}{\@linenumberfalse} + +\newif\if@rotating@loaded \@rotating@loadedfalse +\DeclareOption{rotating}{% + \@rotating@loadedtrue + \AtEndOfClass{\RequirePackage{rotating}}} + +\DeclareOption*{\PassOptionsToClass{\CurrentOption}{article}} +\ProcessOptions* + +\usepackage[T1]{fontenc} +\usepackage{textcomp} +\if@review + \LoadClass[12pt,twoside]{article} + \RequirePackage{setspace} + \doublespacing + \RequirePackage[left=1in,right=1in,top=1in,bottom=1in]{geometry} +\else + \LoadClass[11pt,twoside]{article} + \RequirePackage[left=1.1in,right=1.1in,top=1.1in,bottom=1.1in]{geometry} +\fi + +\if@discussion + \gdef\disfirstpage{0} + \RequirePackage{subfiles} + \if@chinese + \else + \RequirePackage{bibunits} + \defaultbibliographystyle{jds} + \fi +\fi + +\if@chinese + %% By default, scheme = chinese, which changes line space and font size. + %% see the manual of ctex package for details + \RequirePackage[UTF8,heading=true,scheme=plain]{ctex} + \ctexset{ + autoindent = true, + part/afterindent = true, + % chapter/afterindent = true, + section/afterindent = true, + subsection/afterindent = true, + subsubsection/afterindent = true, + paragraph/afterindent = true, + subparagraph/afterindent = true, + contentsname = {??}, + listfigurename = {??}, + indexname = {??}, + abstractname = {??}, + figurename = {?}, + tablename = {?}, + appendixname = {??}, + bibname = {????}, + proofname = {??} + % algorithmname = {??} + } + \RequirePackage{xeCJK} + \RequirePackage{setspace} + \gdef\@pageprefix{S} +\else + \gdef\@pageprefix{} +\fi + + +\def\@@lastpage{\@lastpage} + +\InputIfFileExists{jdsart.cfg}{}{ + \RequirePackage{lastpage} + \def\@@lastpage{\pageref*{LastPage}} + } + +\def\@journal{Journal of Data Science} +\def\@journalabbrev{Journal of Data Science} +\def\@journal@url{https://www.jds.org} +\def\@JOURNAL{JOURNAL OF DATA SCIENCE} +\def\@jid{jds} +\def\@JID{JDS} + +\@ifundefined{doi@base}{\def\doi@base{https://doi.org/}}{} +\def\print@doi{\href{\doi@base\@doi}{\doi@base\@doi}} +\def\@ppages{\@firstpage--\@@lastpage} +\def\print@volume@line{\textsc{\@journal}\space\@volume\space (\@issue), \@ppages\\ \@pubmonth\ \@pubyear} + +\def\@volume{0} +\def\@pubyear{0000} +\def\@pubmonth{} +\def\@@ppages{\empty} +\def\@firstpage{1} +\def\@lastpage{0} +\def\@doi{\empty} + +\def\volume#1{\gdef\@volume{#1}} +\def\issue#1{\gdef\@issue{#1}} +\def\pubyear#1{\gdef\@pubyear{#1}} +\def\pubmonth#1{\gdef\@pubmonth{#1}} +\def\firstpage#1{\gdef\@firstpage{#1}\setcounter{page}{#1}} +\def\lastpage#1{\gdef\@lastpage{#1}} +\def\doi#1{\gdef\@doi{#1}} + +\def\logo@vtex{jds-logo-vtex} +\def\includevtexlogo{\null} +\def\@vtex@url{https://www.vtex.lt/en/} + +% COPYRIGHT +% Define font psyr for copyright: +\DeclareFixedFont{\copyrf}{U}{psy}{m}{n}{8} +\DeclareFixedFont{\copyrfa}{U}{psy}{m}{n}{9} +\def\copyright@sign{{\copyrf\char'323}} + +\newif\if@open@access \global\@open@accesstrue +\def\cclicense@link{https://creativecommons.org/licenses/by/4.0/} +\def\@openaccess{Open access article under the \href{\cclicense@link}{CC BY} license.} +\def\copyright@statement{\textcopyright\ \@pubyear\ The Author(s). Published by the School of Statistics and + the Center for Applied Statistics, Renmin University of China. \@openaccess} +\def\copyright@line{\@copyright@fmt{\copyright@statement}} + +\newbox\copyright@box +\def\@copyright@fmt#1{% + \setbox\copyright@box=\vbox to0pt{\vss \footnotesize #1\hfill}% + \wd\copyright@box=\z@ + \ht\copyright@box=\z@ + \dp\copyright@box=\z@ + \footnotesize + \raise %11pt + \baselineskip + \box\copyright@box} + +\def\nocopyrightline{\gdef\copyright@line{\@copyright@fmt{}}} +\def\credit@line{\@credit@fmt{\@credit}} +\def\@credit@fmt#1{{\parindent0pt#1}} +\def\@credit{% + \vbox to 0pt{% + \vskip-1.85pc + \hskip-\textwidth + \noindent + \raise4mm\hbox to \textwidth{% + \footnotesize + \href{\@journal@url}{www.jds.org}% + \hfill + \href{\@vtex@url}{\includevtexlogo}% + }% + }% + } + +\let\@received\@empty +\let\@revised\@empty +\let\@accepted\@empty +\let\@publishedonline\@empty +\def\received#1{\gdef\@received{{Received on~\hbox{#1}}}} +\def\revised#1{\gdef\@revised{{Revised on~\hbox{#1}}}} +\def\accepted#1{\gdef\@accepted{{Accepted on~\hbox{#1}}}} +\def\publishedonline#1{\gdef\@publishedonline{{Published~online:~\hbox{#1}}}} +\def\print@history{% + \def\sep{\def\sep{, \allowbreak}}% + \ifx\@received\@empty\else\sep\@received\fi% + \ifx\@revised\@empty\else\sep\@revised\fi% + \ifx\@accepted\@empty\else\sep\@accepted\fi% + \ifx\@publishedonline\@empty\else\sep\@publishedonline\fi% + } + +\RequirePackage{amsmath} +\RequirePackage{amssymb} + +\ifdefined\noaccents@ + \noaccents@ +\fi + +\let\kernel@bordmatrix@left( +\let\kernel@bordmatrix@right) +\def\kernel@bordermatrix#1{% + \begingroup + \@nameuse{boredematrix@begin@hook}% + \m@th \@tempdima 8.75\p@ + \setbox\z@\vbox{% + \def\cr{\crcr\noalign{\kern2\p@\global\let\cr\endline}}% + \ialign{$##$\hfil\kern2\p@\kern\@tempdima&\thinspace\hfil$##$\hfil + &&\quad\hfil$##$\hfil\crcr + \omit\strut\hfil\crcr\noalign{\kern-\baselineskip}% + #1\crcr\omit\strut\cr + }% + }% + \setbox\tw@\vbox{\unvcopy\z@\global\setbox\@ne\lastbox}% + \setbox\tw@\hbox{\unhbox\@ne\unskip\global\setbox\@ne\lastbox}% + \setbox\tw@\hbox{$\kern\wd\@ne\kern-\@tempdima + \left\kernel@bordmatrix@left\kern-\wd\@ne + \global\setbox\@ne\vbox{\box\@ne\kern2\p@}% + \vcenter{\kern-\ht\@ne\unvbox\z@\kern-\baselineskip} + \,\right\kernel@bordmatrix@right$% + }% + \null\;\vbox{\kern\ht\@ne\box\tw@}% + \@nameuse{boredematrix@end@hook}% + \endgroup + } +\def\set@bordermatrix@braces#1#2{% + \let\kernel@bordmatrix@left#1\relax + \let\kernel@bordmatrix@right#2\relax + \kernel@bordermatrix + } +\def\bordermatrix#1#{% + \ifx\relax#1\relax + \expandafter\kernel@bordermatrix + \else + \expandafter\set@bordermatrix@braces + \fi#1% + } + +\AtBeginDocument{% + \def\tfrac#1#2{{\textstyle\frac{#1}{#2}}}% + \def\dfrac#1#2{{\displaystyle\frac{#1}{#2}}}% + }% +\def\dvt{\colon\ } +\def\dvtx{\nobreak\,\colon\relax} +% +\let\ge\geqslant +\let\geq\geqslant +\let\le\leqslant +\let\leq\leqslant + +\RequirePackage{amsthm} + +\RequirePackage{graphicx} +\RequirePackage{array} + +%% frontmatter +\newlength\printhistorymargin +\setlength\printhistorymargin{0pc} +\newenvironment{printhistory} + {% + \list{}{% + \rightmargin\printhistorymargin + \leftmargin\printhistorymargin + }% + \item\relax\centering + }{% + \par + \endlist + } + +\newif\iflongmktitle \longmktitlefalse +\def\@@@Label#1{% + \@bsphack + \protected@write\@auxout{}{\string\Newlabel{#1}{\@currentlabel}}% + \@esphack + } +\def\Newlabel#1#2{% + \expandafter\gdef\csname X@#1\endcsname{#2}% + } +\def\@@@Ref#1{% + \@ifundefined{X@#1} + {0} + {\csname X@#1\endcsname}% + } +\RequirePackage{pifont} + +\def\abstract@preskip{28pt} +\newcommand\abstract@fmt[1]{\textbf{#1}} +\newcommand\abstract@name{Abstract} +\newcommand\abstract@name@skip{\vskip7pt} +\newbox\absbox +\renewenvironment{abstract} + {% + \global\setbox\absbox=\vbox\bgroup + \hsize=\textwidth + \def\baselinestretch{1}% + \centerline{\abstract@fmt{\abstract@name}}% + \noindent\unskip + \abstract@name@skip + \noindent\ignorespaces + }{% + \par + \egroup + } + +\newcount\keyboxcnt +\newenvironment{keywords}[1][]% + {% + \expandafter\newsavebox\csname keybox#1\endcsname + \@ifundefined{keybox@list} + {% + \csxdef{keybox@list}{keybox#1}% + }{% + \csxappto{keybox@list}{,keybox#1}% + }% + \def\@tempa{#1}% + \def\PACS@name{PACS}% + \def\JEL@name{JEL}% + \def\MSC@MMX@name{MSC2010}% + \ifx\@tempa\@empty + \def\keywords@prefix{\textbf{Keywords}\enskip\space}% + \else + \ifx\@tempa\PACS@name + \def\keywords@prefix{\it PACS:\space}% + \else + \ifx\@tempa\JEL@name + \def\keywords@prefix{\textbf{JEL~classification}\enskip\space}% + \else + \ifx\@tempa\MSC@MMX@name + \def\keywords@prefix{\textbf{2010~MSC}\enskip\space}% + \def\MSC@active@keywords{1}% + \else + \def\keywords@prefix{\textbf{#1~MSC}\enskip\space}% + \def\msc@tempa{2010}% + \ifx\@tempa\msc@tempa + \def\MSC@active@keywords{1}% + \fi + \fi + \fi + \fi + \fi + \@firstkwdtrue + \global\advance\keyboxcnt\@ne + \global\expandafter\setbox\csname keybox#1\endcsname=\vbox\bgroup + \hsize=\textwidth + \normalfont\normalsize + \def\baselinestretch{1}% + \parskip\z@ + % 10pt between abstract/keyword boxes and baselineskip + 2pt between keyword baselines + \ifdefined\first@keyword@env + \prevdepth\z@ + \medskip + \fi + \noindent{\keywords@prefix}% + \raggedright % Keywords are not justified. + \ignorespaces + }{% + \kwd@lastsep\par + \egroup + \gdef\first@keyword@env{1}% + } + +\ifdefined\directlua + \def\fix@italic@cor{% + %% remove italic correction at the end of inline math before keyword sep ';' + \directlua{ + local prevnode = tex.nest[tex.nest.ptr].tail + if prevnode and prevnode.id == node.id("math") then + local pprevnode = prevnode.prev + if pprevnode and pprevnode.id == node.id("kern") then + node.remove(tex.nest[tex.nest.ptr].head, pprevnode) + end + end + }} +\else + \let\fix@italic@cor\empty +\fi + +\newif\if@firstkwd +\def\kwd@sep{\unskip\fix@italic@cor; } +\def\kwd@lastsep{\unskip} +\def\kwd@font{\itshape} +\def\msc@hrefprefix{https://www.ams.org/msc/msc2010.html?s=} +\newcommand\kwd[2][]{% + \bgroup + \kwd@font + \def\@tempa{#1}% + \if@firstkwd + \global\@firstkwdfalse + \let\kwd@sep\@empty + \fi + \ifx\@tempa\@empty + \kwd@sep + \else + \@tempa + \fi + \ifdefined\MSC@active@keywords + \href{\msc@hrefprefix #2}{#2}% + \else + #2% + \fi + \egroup + \ignorespaces + } + +\newdimen\Columnwidth +\Columnwidth=\columnwidth + +\def\printWarning{% + \mbox{}\par\vfill\par + \bgroup + \fboxsep=12pt + \fboxrule=1pt + \hspace*{.18\textwidth}% + \fcolorbox{gray50}{gray10}{\box\warnbox}% + \egroup + \par\vfill\thispagestyle{empty}% + \setcounter{page}{0}% + \clearpage + } + +\long\def\finalMaketitle{% + % \enlargethispage{-4.5pc}% + \if@discussion + %% for discussion, only show first page header style once + \if\disfirstpage0 + \thispagestyle{copyright} + \gdef\disfirstpage{1} + \xdef\@discussionfirstpage{\thepage} + \else + \if@chinese + \gdef\disfirstpage{2} + \thispagestyle{chfirstpage} + \else + \thispagestyle{doipage} + \fi + \fi + \else + \newpage + \thispagestyle{copyright} + \fi + \MaketitleBox + } + +\def\output@frontmatter@authors{% + \def\do##1{% + \expandafter\author@fmt\expandafter{\the\c@author}{\csuse{##1t@loc@notes}}{##1}% + }% + \front@author@defs + \expandafter\docsvlist\expandafter{\author@list}% + \par + } + +\def\validate@frontmatter@authors{% + \def\do##1{% + \setbox\@tempboxa=\hbox{% + \let\bio\@gobble + \let\degs\@gobble + \let\prefix\@gobble + \let\inits\@gobble + \let\fnms\@gobble + \let\snm\@gobble + \let\suffix\@gobble + \let\roles\@gobble + \def\thanks##1##{\@gobble}% + \def\thanksref##1##{\@gobble}% + \let\thankstext\@gobble + \def\ead####1####{\global\stepcounter{countead}}% + \csname ##1\endcsname}% + }% + \expandafter\docsvlist\expandafter{\author@list}% + } + +\def\output@publishedonline{% + \ifx\@empty\@publishedonline + \else + \publishedonline@prefix\hbox{\@publishedonline}\publishedonline@postfix + \@historytrue + \fi + } + +\AtBeginDocument{\csname \@articletype @articletype@control\endcsname} + +\def\printaug{% + \begin{flushright} + \def\do##1{\csuse{##1}\\}% + \expandafter\docsvlist\expandafter{\author@list}% + \end{flushright} + } + +\ifdefined\HCode + \let\printaug\relax +\else + \def\editorial@articletype@control{% + \let\output@frontmatter@authors\validate@frontmatter@authors + \let\output@publishedonline\relax + } + \def\other@articletype@control{% + \let\output@publishedonline\relax + \gdef\author@list{}% + \nocopyrightline + } +\fi + +\def\address@list{} +\def\title@preskip{0pt} +\def\title@afterskip{15pt} +\def\subtitle@afterskip{18pt} +\let\author@sep\relax +\long\def\MaketitleBox{% + \def\baselinestretch{1}% + \begin{center}% + \def\baselinestretch{1}% + %\null + \vskip\title@preskip + \Large\textbf{\title@fmt{\@title}{\csuse{title@t@loc@notes}}}% + \par + \vskip\title@afterskip + \ifdefined\@subtitle + \large\textbf{\subtitle@fmt{\@subtitle}{\csuse{subtitle@t@loc@notes}}}% + \par + \vskip\subtitle@afterskip + \fi + \output@frontmatter@authors + \par + \def\do##1{% + \address@fmt{\csuse{##1@counter}}{\csuse{##1t@loc@notes}}{\csuse{##1}}% + }% + \expandafter\docsvlist\expandafter{\address@list}% + \end{center}% + \ifdefined\@dedicated + \unskip\vskip9.5pt + \dedicated@fmt{\@dedicated}% + \par\vskip15pt + \fi + \ifvoid\absbox + \else + \unskip + \vskip\abstract@preskip + \unvbox\absbox + \par + \vskip10pt + \fi + \ifdefined\keybox@list + \newdimen\saveboxdepth + \count@=\z@ + \def\do##1{% + \saveboxdepth=\expandafter\dp\csname##1\endcsname + \expandafter\unvbox\expandafter\csname ##1\endcsname + \advance\count@\@ne + % depth is 0 for all keyword boxes except the last one + \ifnum\count@=\keyboxcnt + \prevdepth=\the\saveboxdepth + \else + \addvspace{-\the\saveboxdepth}% + \fi + }% + \expandafter\docsvlist\expandafter{\keybox@list}% + \par + \vskip10pt + \fi + } + +\def\FNtext#1{\par\bgroup\footnotesize#1\egroup} +\newdimen\space@left +\def\alarm#1{% + \typeout{******************************}% + \typeout{#1}% + \typeout{******************************}% + } +\long\def\getSpaceLeft{% + \global\setbox0=\vbox{\hsize=\textwidth\MaketitleBox}% + \global\setbox1=\vbox{\hsize=\textwidth + \let\footnotetext\FNtext + \printFirstPageNotes}% + \xdef\noteheight{\the\ht1}% + \xdef\titleheight{\the\ht0}% + \@tempdima=\vsize + \advance\@tempdima-\noteheight + \advance\@tempdima-1\baselineskip + } +\skip\footins=24pt +\newbox\@@@@boxa +\newbox\@@@@boxb +\def\maketitle{% + \iflongmktitle\getSpaceLeft + \global\setbox\@@@@boxa=\vsplit0 to \@tempdima + \box\@@@@boxa\par\resetTitleCounters + \thispagestyle{copyright}% + \printFirstPageNotes + \box0% + \else + \finalMaketitle + \fi + }% + +\def\runtitle#1{\gdef\@runtitle{\hbox{#1}}} +\def\@runtitle{} +\def\runauthor#1{\gdef\@runauthor{#1}} +\def\@runauthor{\output@head@authors} +\def\output@head@authors{} + +\def\runauthor@etal{~et al.} + +\def\head@author@defs{% + \front@author@defs + \let\write@aux@au\@gobbletwo + \def\ead##1##{\@gobble}% + \def\@inits##1{}% + \def\@snm##1{}% + \def\@fnm##1{}% + \let\~\empty + \def\inits##1{\def\@inits{, ##1}}% + \def\snm##1{\def\@snm{##1}}% + \def\fnms##1{\def\@fnms{##1}\ignorespaces}% + } + +\def\head@author@fmt#1{% + \head@author@defs + \get@author@num #1\relax + \set@author@sep + \@author@sep + \@nameuse{#1}% + \@snm \@inits + } +\def\head@authors@font{\itshape\selectfont} + +\newbox\headbox +\def\@output@head@authors{% + \bgroup + \head@authors@font + \ifnum\c@author>2\relax + \head@author@fmt{author1}\runauthor@etal + \else + \def\do##1{\head@author@fmt{##1}}% + \expandafter\docsvlist\expandafter{\author@list}% + \fi + \egroup + } + +%% define header and footer +\RequirePackage{fancyhdr} +\setlength{\headheight}{24pt} +\pagestyle{fancy} +\fancyhead{} % clear all header fields +\fancyhead[CO]{% + \csname check@head@line\endcsname% + \textit \@runtitle% +} +\fancyhead[CE]{\@runauthor} +\fancyhead[LE,RO]{\@pageprefix\thepage} +\fancyfoot{} % clear all footer fields +\renewcommand{\headrulewidth}{0pt} +\renewcommand{\footrulewidth}{0pt} + +\def\check@head@line{% + \def\pseudo@head@line{\@pageprefix 000\space{\textit\@runtitle}\space\@pageprefix 000}% + \bgroup + \settowidth\@tempdima{\pseudo@head@line}% + \ifdim\@tempdima>\textwidth\relax + \ClassError{jdsart}{Running title is too long!}{}% + \global\let\check@head@line\relax + \fi + \egroup +} + +\fancypagestyle{copyright}{% + \renewcommand{\headrulewidth}{0pt}% + \fancyhf{}% + \fancyhead[L]{\normalfont\footnotesize\print@volume@line} + \fancyhead[R]{\footnotesize DOI: \@doi\\ \@pretitle} + \fancyfoot[LE,RO]{}% number of the first page + \fancyfoot[RE,LO]{\footnotesize \history@fmt} +} + +\if@discussion + \fancypagestyle{doipage}{% + \renewcommand{\headrulewidth}{0pt}% + \fancyhf{}% + \fancyhead[L]{\footnotesize {\scshape \@journal} + \@volume\ (\@issue), + \IfRefUndefinedExpandable{\@doi}{ + \@pageprefix\thepage--\pageref*{LastPage}\\ + }{ + \@pageprefix\thepage--\pageref*{\@doi}\\ + } + \@smonth\ \@syear + } + \fancyhead[R]{\footnotesize DOI:\ \@doi\\ \@pretitle} + \fancyfoot[LE,RO]{\@pageprefix\thepage}% number of the page + } +\fi + +\if@chinese + \fancypagestyle{chfirstpage}{% + \renewcommand{\headrulewidth}{0pt}% + \fancyhf{}% + \fancyhead[L]{\footnotesize {\scshape \@journal} + \@volume\ (\@issue), + \@pageprefix\@discussionfirstpage--\@pageprefix\pageref*{LastPage}\\ + \@syear 年 \@smonth 月 + } + \fancyhead[R]{\footnotesize DOI:\ \@doi\\ \@pretitle} + \fancyfoot[LE,RO]{\@pageprefix\thepage}% number of the first page + \fancyfoot[RE,LO]{\footnotesize \history@fmt} + } +\fi + +\if@linenumber + \long\gdef\numberlines@size{\scriptsize\ttfamily} + \long\gdef\numberlines@skip{24\p@} + % + \def\numberlines@hook{% + \appto\@evenhead\put@numberlines@box + \appto\@oddhead\put@numberlines@box + } + % + \gappto\ps@copyright\numberlines@hook + \numberlines@hook + % + \newbox\numberlines@box + \newskip\numberlines@box@skip + % + \def\set@numberlines@box{% + \setlength\numberlines@box@skip\headsep + \addtolength\numberlines@box@skip{5\p@}% + % + \setbox\numberlines@box\vtop to\textheight{% + \parindent=\z@ + \vskip\z@ + \@tempcnta=0% + \@tempdima=\z@ + \loop + \advance\@tempcnta by 1% + \advance\@tempdima by \baselineskip + \hbox to\textwidth{% + \llap{\numberlines@size\the\@tempcnta\kern\numberlines@skip}% + \hfill + \rlap{\numberlines@size\kern\numberlines@skip\the\@tempcnta}% + }% + \ifdim\@tempdima<\textheight\repeat + \vss + }% + % + \ht\numberlines@box=\z@ + \dp\numberlines@box=\z@ + } + % + \def\put@numberlines@box{\lower\numberlines@box@skip\hbox to\z@{\hss\copy\numberlines@box}} + % + \AtBeginDocument{\set@numberlines@box} + % +\fi + +%% caption +% \RequirePackage{caption} +% \captionsetup{margin=0pt,font=small,labelfont=bf,labelsep=period} +% \captionsetup[table]{position=top} +% \renewcommand\figurename{Fig.} + +\def\floatfmt{\centering} +\appto\@floatboxreset\floatfmt + +\newtoks\fm@temp@toks +\newenvironment{frontmatter} + {% + \global\@firstauthortrue + \let\saved@thefootnote\thefootnote + }{% + \do@pdfmeta@authors + \maketitle + \csname check@articletype@pretitle@status\endcsname + } + +\def\init@standard@footnote{% + \ifdefined\init@standard@footnote@done + \else + \gdef\init@standard@footnote@done{1}% + \global\let\thefootnote\saved@thefootnote + \setcounter{footnote}{0}% + \fi + } + +%% lists +\setlength{\parskip}{0ex plus0.1ex minus0.1ex} +\setlength{\parindent}{2em} + +\RequirePackage{enumitem} +\setlist{nosep, leftmargin=1.5em, labelwidth=0.25em} +\raggedbottom + +\IfFileExists{vtexnatbib.sty} + {\RequirePackage[authoryear,round]{vtexnatbib}} + {\RequirePackage[authoryear,round]{natbib}} + +\let\code\texttt +\let\proglang\textsf +\newcommand{\pkg}[1]{{\fontshape{it}\selectfont #1}} + +\setlength{\bibsep}{0em} + +\def\MR@url{https://www.ams.org/mathscinet-getitem?mr=} +\def\MR#1{\href{\MR@url#1}{MR#1}} + +\RequirePackage{xcolor} +\definecolor{linkblue}{rgb}{0 0.4235 0.6392} +\RequirePackage{hyperref} +\hypersetup{% + colorlinks=true, + pdfborder={0 0 1}, + linkcolor=linkblue, + citecolor=linkblue, + urlcolor=linkblue, + breaklinks, + hypertexnames=false + } + + + +%% >>> startlocaldefs, endlocaldefs >>> +\let\check@localdefs@status\relax +\def\startlocaldefs{% + \gdef\check@localdefs@status{\dbl@error{Missing command \string\endlocaldefs!}}% + \gdef\startlocaldefs{\dbl@error{Command \string\startlocaldefs already exists!}}% + \makeatletter} +\def\endlocaldefs{% + \makeatother + \ifx\check@localdefs@status\relax + \gdef\check@localdefs@status{\dbl@error{Missing command \string\startlocaldefs!}}% + \else + \global\let\check@localdefs@status\relax + \fi + \gdef\endlocaldefs{\dbl@error{Command \string\endlocaldefs already exists!}}% + } +\AtBeginDocument{\check@localdefs@status}% +%% <<< startlocaldefs, endlocaldefs <<< + +%% >>> \bid >>> +\define@key{bid}{mr}{\def\bid@mr{\MR{#1}}} +\define@key{bid}{doi}{\def\bid@doi{\doiurl{#1}}} +\define@key{bid}{pubmed}{} +\define@key{bid}{pii}{} + +\def\bid#1{% + \bgroup + \let\bid@output\@empty + \setkeys{bid}{#1}\ignorespaces% + % Output sequence: + % . + \ifdefvoid\bid@mr{}{\appto\bid@output{\bid@mr}}% + \ifdefvoid\bid@doi{}{% + \ifdefempty\bid@output{}{\appto\bid@output{. }}% separator if output not empty + \appto\bid@output{\bid@doi}% + }% + \bid@output + \egroup + } + +\def\set@proof@mr{% + \define@key{bid}{mr}{% + \@nameuse{startchkuriskip}% + \safe@href{https://www.ams.org/mathscinet-getitem?mr=##1}{MR##1}% + \@nameuse{endchkuriskip}% + \ignorespaces + }% + \define@key{bid}{doi}{% + \ifnum\doc@stage<200\relax + \@nameuse{startchkuriskip}% + \safe@href{\doi@base ##1}{\doi@base ##1}% + \@nameuse{endchkuriskip}% + \space + \fi + }% + } +\def\set@crc@mr{% + \define@key{bid}{mr}{% + \if@sici + \global\@biblinktrue + \fi + \ignorespaces + }% + } +%% <<< \bid <<< + +\frenchspacing + +%% >>> Metadata >>> +\def\@publisher{School of Statistics, Renmin University of China} +\def\@publisheraddress{No. 59 Zhongguancun Street, Haidian + District Beijing 100872, P.R. China} +\def\@issn{1680-743X} +\def\@eissn{1683-8602} +\def\@articletype{} +\def\articletype#1{\gdef\@articletype{#1}} +\@onlypreamble\articletype +\def\aid#1{\gdef\@aid{#1}} +%% <<< Metadata <<< + +%% >>> Frontmatter >>> +\def\@pretitle{} +\def\pretitle#1{\gdef\@pretitle{#1}} +\newenvironment{aug} + {} + {% + \global\let\thefootnote\saved@thefootnote + \setcounter{footnote}{\the\c@authornote}% footnote counter is continuous with authornote + } + +\ifdefined\HCode +\else + \def\lmr@font@fix{% + \DeclareFontShape{T1}{lmr}{b}{sc}{<->ssub*cmr/bx/sc}{} + \DeclareFontShape{T1}{lmr}{bx}{sc}{<->ssub*cmr/bx/sc}{} + } + \@nameuse{\rmdefault @font@fix} +\fi + +%% >>> author setup >>> +\def\author@font{\fontsize{\@xipt}{17pt}\scshape} +\def\author@preskip{2\p@} +\newif\if@firstauthor +\newif\if@newelem +\let\cmd@for@aut@elem\@gobblethree +\let\mdforvkt\@gobblethree +\let\orcid\@gobble + +\newcounter{author} + +\protected\def\fmt@snm#1{\textsc{#1}} +\protected\def\fmt@fnms#1{\textsc{#1}} + +\let\write@aux@au\@gobbletwo +\protected\def\@write@aux@au#1#2{% + \protect\protected@write\@auxout{}{\string\global\string\@namedef{metadata@\author@num @#1}{#2}}% + } + +\def\front@author@eval@emails{% + \let\bio\@gobble + \let\fnms\@gobble + \let\snm\@gobble + \let\inits\@gobble + \let\prefix\@gobble + \let\degs\@gobble + \let\suffix\@gobble + \let\roles\@gobble + \let\~\empty + } + +\def\front@author@defs{% + \let\bio\@gobble + \let\write@aux@au\@write@aux@au + \protected\def\degs##1{##1}% + \protected\def\prefix##1{% + \write@aux@au{prefix}{##1}% + \fmt@prefix{##1}% + }% + \protected\def\fnms##1{\fmt@fnms{##1}}% + \protected\def\inits##1{% + \def\@inits{##1}% + \write@aux@au{inits}{##1}% + }% + \protected\def\snm##1{% + \def\@snm{##1}% + \write@aux@au{snm}{##1}% + \fmt@snm{##1}% + }% + \protected\def\suffix##1{##1}% + \protected\def\roles##1{##1}% + } + +\def\author{% + \@ifnextchar[% + {\author@optarg} + {\author@optarg[]}% + } + +\def\author@optarg[#1]#2{% + %% headings will use authors + \global\let\output@head@authors\@output@head@authors + \stepcounter{author}% + \beg@elem + \@for\@tempa:=#1\do{\expandafter\add@addressref\expandafter{\@tempa}}% + \let\oldunderscore\_% + \let\_\textunderscore + \init@fm@element{author}{#2}% + \let\_\oldunderscore + \csxdef{author\the\c@author}{\expandonce{\@author}}% + \protected\csxdef{author\the\c@author t@loc@notes}{\the\t@loc@notes}% + \@ifundefined{author@list} + {% + \xdef\author@list{author\the\c@author}% + }{% + \xappto\author@list{,author\the\c@author}% + }% + \bgroup + \front@author@eval@emails + \setbox\z@=\hbox{#2}% + \egroup + \csname author@end@hook\endcsname + }% + +\newtoks\pdfmeta@authors +\def\addto@pdfmeta@authors#1{% + \begingroup + \xdef\@act{\global\noexpand\pdfmeta@authors{\the\pdfmeta@authors#1}}% + \@act + \endgroup + } + +\def\do@pdfmeta@authors{% + \@tempcnta\z@ + \@whilenum\@tempcnta<\c@author\do{% + \advance\@tempcnta\@ne + \ifnum\@tempcnta>\@ne + \addto@pdfmeta@authors{, }% + \fi + \addto@pdfmeta@authors{% + \@ifundefined{metadata@\the\@tempcnta @inits}{}{% + \@nameuse{metadata@\the\@tempcnta @inits}% + \space + }% + \@nameuse{metadata@\the\@tempcnta @snm}% + }% + }% + \hypersetup{% + pdfinfo={% + Title={\@title}, + Subject={}, + Author={\the\pdfmeta@authors}, + Creator={VTeX}, + Producer={VTeX}, + }% + }% + } + +\def\author@num{0} +\def\get@author@num author#1\relax{\def\author@num{#1}} +\def\author@sep{, } +\def\author@twosep{ and } +\def\author@lastsep{, and } + +\def\set@author@sep{% + \ifnum \author@num>1\relax + \ifnum \author@num=\the\c@author\relax + \ifnum\author@num=2\relax + \let\@author@sep\author@twosep + \else + \let\@author@sep\author@lastsep + \fi + \else + \let\@author@sep\author@sep + \fi + \else + \let\@author@sep\empty + \fi + } + +\newtoks\temptoks +\def\author@fmt#1#2#3{% + \ifdefined\author@fmt@init@done + \else + \author@fmt@init + \fi + \bgroup + \author@font + \get@author@num #3\relax + \set@author@sep + \@author@sep + \@nameuse{#3}% + \edef\@tempb{#2}% + \ifx\@tempb\@empty + \else + \fm@textsuperscript{#2}% + \fi + \egroup + } +\def\fm@textsuperscript#1{% + %$^{\mathrm{\scriptsize #1}}$% + \unskip\textsuperscript{{\normalfont #1}}% + } + +\def\author@fmt@init{% + \gdef\author@fmt@init@done{1}% + \vskip\author@preskip + \leavevmode + } + +\def\author@note@fmt{% + \def\thefootnote{\sxarabic{footnote}}} + +\def\sxarabic#1{% + \expandafter + \ifcase\value{#1}% + \or *% + \or **% + \or ***% + \or ****% + \or *****% + \fi + } + +\def\add@addressref#1{% + \global\advance\note@cnt\@ne + \ifnum\note@cnt>\@ne + \add@xtok\t@loc@notes{\note@sep}% + \fi + \add@tok\t@loc@notes{\ref{aff.#1}}% + } +\def\note@sep{,} +%% <<< author setup <<< + +%% >>> affil setup >>> +\newcounter{address} +\def\address@font{\mdseries\slshape\small} +\def\address@preskip{5\p@} +\def\address{% + \@ifstar + {\address@star} + {\@ifnextchar[{\address@optarg}{\address@noptarg}}% + } + +\def\address@optarg[#1]#2{% + \refstepcounter{address}% + \beg@elem + \init@fm@element{address}{#2}% + \csxdef{address\the\c@address @counter}{\the\c@address}% + \temptoks=\expandafter{\@address}% + \csxdef{address\the\c@address}{\the\temptoks}% + \protected\csxdef{address\the\c@address t@loc@notes}{\the\t@loc@notes}% + \@ifundefined{address@list} + {% + \xdef\address@list{address\the\c@address}% + }{% + \xappto\address@list{,address\the\c@address}% + }% + \label{aff.#1}% + \csname address@end@hook\endcsname + \ignorespaces + } + +\def\address@noptarg#1{% + \refstepcounter{address}% + \beg@elem + \init@fm@element{address}{#1}% + \csxdef{address\the\c@address @counter}{\z@}% + \temptoks=\expandafter{\@address}% + \csxdef{address\the\c@address}{\the\temptoks}% + \protected\csxdef{address\the\c@address t@loc@notes}{\the\t@loc@notes}% + \@ifundefined{address@list} + {% + \xdef\address@list{address\the\c@address}% + }{% + \xappto\address@list{,address\the\c@address}% + }% + \csname address@end@hook\endcsname + \ignorespaces + } + +\def\address@star#1{% + \beg@elem + \init@fm@element{address}{#1}% + \address@fmt{\m@ne}{\the\t@loc@notes}{\@address}% + \ignorespaces + } + +\def\theaddress{\the\c@address} + +\def\address@fmt#1#2#3{% + \ifdefined\address@fmt@init@done + \else + \address@fmt@init + \fi + \bgroup + \def\baselinestretch{1}% + \address@font + \ifnum#1=\z@ + #3\,\fm@textsuperscript{#2}\space + \else + \ifnum#1=\m@ne + \phantom{\fm@textsuperscript{#1}\,}#3\,\fm@textsuperscript{#2}% + \else + \fm@textsuperscript{#1}#3\fm@textsuperscript{#2}% + \fi + \fi + \par + \egroup + } +\def\address@fmt@init{% + \gdef\address@fmt@init@done{1}% + \vskip\address@preskip + \leavevmode + } +%% <<< affil setup <<< + +%% >>> \beg@elem, \proc@elem (based on iosart2x.cls) >>> +\def\beg@elem{% + \global\t@loc@notes={}% + \global\note@cnt\z@ + } +\def\@xnamedef#1{% + \expandafter\xdef\csname #1\endcsname + } +\def\no@harm{% + \let\\=\relax \let\rm\relax + \let\ss=\relax \let\ae=\relax \let\oe=\relax + \let\AE=\relax \let\OE=\relax + \let\o=\relax \let\O=\relax + \let\i=\relax \let\j=\relax + \let\aa=\relax \let\AA=\relax + \let\l=\relax \let\L=\relax + \let\d=\relax \let\b=\relax \let\c=\relax + \let\bar=\relax + \def\qq@group@start{}\def\qq@group@end{}\let\qq@cmd\@firstofone\def\qq@style{}% + \def\protect{\noexpand\protect\noexpand}% + } + +\def\thanks#1#{\@gobble} +\def\thanksref#1#{\@gobble} + +\def\proc@elem#1#2{% + \begingroup + \no@harm + \def\thanks##1##{\@gobble}% + \def\thanksref##1##{\@gobble}% + \csgdef{@#1}{#2}% + \endgroup + \prev@elem=\cur@elem + \cur@elem=\csname e@#1\endcsname + \elem@nothanks#2\thanks\relax + \elem@nothanksref#2\thanksref\relax + } +%% <<< \beg@elem, \proc@elem (based on iosart2x.cls) <<< + +%% >>> thanks macro block >>> +\newtoks\t@glob@notes +\newtoks\t@loc@notes +\newcount\note@cnt +\newcount\c@ThanksID +\newcount\prev@elem \prev@elem=0 +\newcount\cur@elem \cur@elem=0 +\chardef\e@pretitle=1 +\chardef\e@title=1 +\chardef\e@subtitle=1 +\chardef\e@author=2 +\chardef\e@address=3 + +\newif\if@thankshyperlink \global\@thankshyperlinkfalse + +\def\add@tok#1#2{% + \global#1\expandafter{\the#1#2}% + } +\def\add@xtok#1#2{% + \begingroup + \no@harm + \xdef\@act{% + \global\noexpand#1{\the#1#2}% + }% + \@act + \endgroup + } + +\def\@fnsymbol#1{% + \ifcase #1% + \or \TextOrMath \textdagger \dagger + \or \TextOrMath \textdaggerdbl \ddagger + \or \TextOrMath \textsection \mathsection + \or \TextOrMath \textparagraph \mathparagraph + \or \TextOrMath \textbardbl \|% + \or \TextOrMath {\textasteriskcentered \textasteriskcentered }{**}% + \or \TextOrMath {\textdagger \textdagger }{\dagger \dagger }% + \or \TextOrMath {\textdaggerdbl \textdaggerdbl }{\ddagger \ddagger }% + \else \@ctrerr \fi} + +\newcounter{correspnote} +\def\corresp@note@counter{correspnote} + +\newcounter{titlenote} +\def\title@note@counter{titlenote} +\def\subtitle@note@counter{titlenote} +\def\relatedarticle@note@counter{titlenote} +\def\correctionnote@note@counter{titlenote} +\def\correctedarticle@note@counter{titlenote} +\def\commentaryarticle@note@counter{titlenote} + +\newcounter{authornote} +\def\author@note@counter{authornote} +\def\theauthornote{\fnsymbol{authornote}} +\def\thecorrespnote{% + \ifcase\c@correspnote\relax + \or + $\ast$% + \or + $\ast\ast$% + \or + $\ast\ast\ast$% + \fi + }% +\def\thetitlenote{% + \ifcase\c@titlenote\relax + \or + \protect\ding{73}% + \or + \protect\ding{73}\protect\ding{73}% + \or + \protect\ding{73}\protect\ding{73}\protect\ding{73}% + \fi + }% +\let\thesubtitlenote\thetitlenote + +\let\thanks@id\@empty +\let\thanks@mark\@empty +\let\thanks@type\@empty +\define@key{thanks}{id}{\def\thanks@id{#1}} +\define@key{thanks}{mark}{\def\thanks@mark{#1}} +\define@key{thanks}{type}{\def\thanks@type{#1}} +\def\reset@thanks{% + \let\thanks@id\@empty + \let\thanks@mark\@empty + \let\thanks@type\@empty + } +\def\thanks@optarg@fmt#1#2{% + \freefootnotetext[#1]{#2}% + } +\def\thanks@optarg[#1]#2{% + \let\thanks@type\@empty + \setkeys{thanks}{#1}% + \ifx\thanks@type\@empty + \let\thanks@type\fm@object + \fi + \edef\current@counter{\csname\thanks@type @note@counter\endcsname}% + \refstepcounter{\current@counter}% + \setcounter{footnote}{\the\csname c@\current@counter\endcsname}% + \def\thefootnote{\csname the\current@counter\endcsname}% + \edef\thanks@optarg@id{\thanks@type.thanks\the\c@footnote.\jobname}% + \@@@Label{\thanks@optarg@id}% + \if@thankshyperlink %% TODO: make active links + \ifx\thanks@id\@empty% + \Hy@GlobalStepCount\Hy@linkcounter + \xdef\thanks@currentHref{\thanks@type.thanks\the\c@footnote.\jobname.\the\Hy@linkcounter}% + \edef\my@hy@tmp{\noexpand\hypertarget{\thanks@currentHref}{}}% + %\edef\my@hy@tmp{\noexpand\hypertarget{\@currentHref}{}}% in case \refstepcounter{footnote} is used + \else + \edef\my@hy@tmp{\noexpand\hypertarget{\thanks@id thanks}{}}% + \fi + %\expandafter\add@tok\expandafter\t@glob@notes\expandafter{\expandafter{\my@hy@tmp #2}}% + \ifx\thanks@mark\@empty + \def\thanks@mark{\thefootnote}% + \fi + \freefootnotetext[\thanks@mark]{{\my@hy@tmp #2}}% + \else + \ifx\thanks@mark\@empty + \def\thanks@mark{\thefootnote}% + \fi + \thanks@optarg@fmt{\thanks@mark}{#2}% + \fi + \csname end@thanks@optarg@hook\endcsname + \ignorespaces + }% +\def\freefootnotetext[#1]{% + \begingroup + \unrestored@protected@xdef\@thefnmark{#1}% + \endgroup + \@footnotetext + } +%% TODO: do we need \tailthanksref? +\def\tailthanksref[#1]#2{% + \noexpand\pthanksref{#1}% + } +\def\pthanksref#1{% + \global\advance\note@cnt\@ne + \ifnum\note@cnt>\@ne + \global\t@loc@notes\expandafter{\the\t@loc@notes\note@sep}% + \fi + \global\t@loc@notes\expandafter{\the\t@loc@notes#1}% + } + +\def\elem@nothanksref#1\thanksref{% + \futurelet\@peektok\elem@thanksref + } +\def\elem@thanksref{% + \ifx\@peektok\relax + \else + \expandafter\elem@morethanksref + \fi + } +\def\elem@morethanksref#1{% + \csxdef{#1@thankstext@type}{\fm@object}% + \csletcs{thankstext@type}{#1@thankstext@type}% + \add@thanks{id=#1}% + \elem@nothanksref + } + +\def\elem@nothanks#1\thanks{% + \futurelet\@peektok\elem@thanks + } +\def\elem@thanks{% + \ifx\@peektok\relax + \else + \ifx\@peektok[% + \expandafter\expandafter\expandafter\elem@morethankse + \else + \expandafter\expandafter\expandafter\elem@morethanks + \fi + \fi + } +\def\elem@morethankse[#1]#2{% + \thanks@optarg[#1]{#2}% + \add@thanks{#1}% + \elem@nothanks + } +\def\elem@morethanks#1{% + \thanks@optarg[]{#1}% + \add@thanks{}% + \elem@nothanks + } + +\define@key{thankstext}{type}{\gdef\thankstext@type{#1}} +\define@key{thankstext}{id}{\gdef\thankstext@id{#1}} +\def\reset@thankstext{% + \let\thankstext@id\@empty + \let\thankstext@mark\@empty + \let\thankstext@type\@empty + } +\def\thankstext@fmt#1#2{% + \freefootnotetext[#1]{#2}% + } +\def\thankstext[#1]#2{% + \reset@thankstext + \setkeys{thankstext}{#1}% + \ifx\thankstext@type\@empty + \csletcs{thankstext@type}{\thankstext@id @thankstext@type}% + \fi + \edef\current@counter{\csname\thankstext@type @note@counter\endcsname}% + \refstepcounter{\current@counter}% + \setcounter{footnote}{\the\csname c@\current@counter\endcsname}% + \def\thefootnote{\csname the\current@counter\endcsname}% + \@@@Label{\thankstext@id}% + \thankstext@fmt{\thefootnote}{\csname \thankstext@type @prehook\endcsname #2\csname \thankstext@type @posthook\endcsname}% + } + + +\define@key{addthanks}{id}{\def\addthanks@id{#1}} +\define@key{addthanks}{mark}{\def\addthanks@mark{#1}} +\define@key{addthanks}{type}{\def\addthanks@type{#1}} +\def\reset@addthanks{% + \let\addthanks@id\@empty + \let\addthanks@mark\@empty + \let\addthanks@type\@empty + } +\def\add@thanks#1{% + \reset@addthanks + \setkeys{addthanks}{#1}% + \global\advance\note@cnt\@ne + \ifnum\note@cnt>\@ne + \add@xtok\t@loc@notes{\note@sep}% + \fi + \if@thankshyperlink %% TODO: make active links + \ifx\addthanks@id\@empty + \add@xtok\t@loc@notes{\protect\Hy@breaklinkstrue + \protect\hyper@linkstart{link}{\thanks@currentHref}\thefootnote\protect\hyper@linkend}% + %\protect\hyper@linkstart{link}{\@currentHref}\thefootnote\protect\hyper@linkend}% + \else + \add@xtok\t@loc@notes{\protect\Hy@breaklinkstrue + \protect\hyper@linkstart{link}{\addthanks@id thanks}\addthanks@id\protect\hyper@linkend}% + \fi% + \else + \ifx\addthanks@id\@empty + \def\addthanks@mark{\thefootnote}% + \else + \ifx\addthanks@mark\@empty + \def\addthanks@mark{\@@@Ref{\addthanks@id}}% + \fi + \fi + \add@xtok\t@loc@notes{\addthanks@mark}% + \fi + \csname end@add@thanks@hook\endcsname + } +\newcount\c@sv@footnote +\def\output@glob@notes{% + \bgroup + \the\t@glob@notes + \egroup + } +%% <<< thanks macro block <<< + +%% >>> init FM element setup >>> +\newcommand\init@fm@element[3][]{% + \def\fm@object{#2}% + \ifcsname init@fm@#2@done\endcsname + \else + \setkeys{#2}{#1}% + \fi + \csgdef{init@fm@#2@done}{1}% + \proc@elem{#2}{#3}% + } +%% <<< init FM element setup <<< + +%% >>> title/subtitle setup >>> + +\def\title#1{% + \beg@elem + \gdef\@runtitle{\hbox{#1}}% + \init@fm@element{title}{#1}% + \protected\csxdef{title@t@loc@notes}{\the\t@loc@notes}% + \csname title@end@hook\endcsname + } +\def\subtitle#1{% + \beg@elem + \init@fm@element{subtitle}{#1}% + \protected\csxdef{subtitle@t@loc@notes}{\the\t@loc@notes}% + \csname subtitle@end@hook\endcsname + } +\def\title@fmt#1#2{% + \bgroup + #1\hbox{\title@textsuperscript{#2}}% + \egroup + } +\let\subtitle@fmt\title@fmt +\def\title@textsuperscript#1{\textsuperscript{#1}} +%% <<< title/subtitle setup <<< + +%% >>> email >>> +\def\ead@mark{} +\def\ead@href#1#2{\href{#1}{\ead@mark}} +\def\email@prefix{Email:\space} +\def\email@postfix{.} +\def\email@font{ } + +\newtoks\email@author@toks +\def\email@author@fmt{\ignorespaces} + +%% blocked use of printead.sty in author environment +\newcounter{countead} +\def\ead@hook#1#2{% + \stepcounter{countead}% + \ifx.#2. + \gdef\ead@mark{???}% + \else + \gdef\ead@mark{#2}% + \fi + } + +\def\ead@hrefprefix@email{mailto:} + +\DeclareRobustCommand\ead[2][]{{% +%% few necessary printead.sty definitions + \define@key{ead}{label}{\def\ead@id{#1}}% + \define@key{ead}{email}[true]{\def\ead@type{email}}% + \def\@setead##1##2##3{% + \bgroup + \csname @setead@hook\endcsname + \expandafter\xdef\csname ead@##1@##2\endcsname{##3}% + \egroup + }% + % + \csname before@ead@hook\endcsname + \setkeys{ead}{label={},email,#1}% + \@setead{\ead@id}{options}{#1}% + \@setead{\ead@id}{type}{\ead@type}% + \@setead{\ead@id}{value}{#2}% + \csedef{\ead@id @email@author}{}% + \def\email@type{email}% + \ifx\ead@type\email@type + \ead@hook{#1}{#2}% + \add@xtok\email@author@toks{% + \csuse{emailauthor@sep}% + \ead@href{\ead@hrefprefix@email #2}{#2}% + }% + \gdef\emailauthor@sep{\unskip\space or\space}% + \fi + \csname after@ead@hook\endcsname + }\ignorespaces} + +\def\output@frontmatter@emails{% + \bgroup + \edef\emaillist{\the\email@author@toks}% + \ifx\emaillist\empty + \else + \bgroup + \email@font + \email@prefix \the\email@author@toks \email@postfix + \egroup + \fi + \egroup + } + +\protected\def\corresp@posthook{\output@frontmatter@emails} +%% <<< email <<< + +%% >>> history >>> +%% history month - full name +\def\setvaluelist#1#2{% + \@tempcnta=0\relax + \@for\@curr@val:=#2\do{% + \advance\@tempcnta by1\relax + \expandafter\protected@xdef\csname #1@item@\the\@tempcnta\endcsname{\@curr@val}% + }% + \expandafter\protected@xdef\csname #1@item@0\endcsname{\the\@tempcnta}% + } +\xdef\getitemvalue#1#2{\noexpand\csname #1@item@#2\endcsname} +\def\sday#1{\gdef\@sday{#1}#1} +\def\syear#1{\gdef\@syear{#1}#1} +\def\smonth#1{% + \gdef\@smonth{#1}% + \@ifundefined{month@item@#1}% + {\@latex@error{There is no month with number '#1'!}\@ehc}% + {\getitemvalue{month}{#1}}% + }% +\setvaluelist{month}{January,February,March,April,May,June,July,August,September,October,November,December} +\expandafter\def\csname month@item@0\endcsname{????} + +\newif\if@history \@historyfalse +\def\define@history@date#1{% + \expandafter\def\csname #1\endcsname##1{% + \def\@tempa{##1}% + \ifx\@tempa\@empty + \else + \csgdef{@#1}{##1}% + \fi + }% + } +\define@history@date{received} +\define@history@date{revised} +\define@history@date{accepted} +\define@history@date{publishedonline} + + +\def\received@prefix{Received~} +\def\revised@prefix{Revised~} +\def\accepted@prefix{Accepted~} +\def\publishedonline@prefix{Published~online:~} +\if@chinese + \def\received@prefix{投稿于~} + \def\accepted@prefix{接收于~} +\fi + +\def\received@postfix{} +\def\revised@postfix{} +\def\accepted@postfix{} +\def\publishedonline@postfix{} +\def\history@prefix{} +\def\history@postfix{} + +\def\history@sep{% + \if@history + ; \allowbreak + \fi + } +\def\history@fmt{% + \history@prefix + \ifx\@empty\@received + \else + \received@prefix\hbox{\@received}\received@postfix + \@historytrue + \fi + \ifx\@empty\@revised + \else + \history@sep + \revised@prefix\hbox{\@revised}\revised@postfix + \@historytrue + \fi + \ifx\@empty\@accepted + \else + \history@sep + \accepted@prefix\hbox{\@accepted}\accepted@postfix + \@historytrue + \fi + \output@publishedonline + \history@postfix + } + +\def\history@skip{14pt} +\def\output@frontmatter@history{% + \vskip\history@skip + \bgroup + \small + \history@fmt + \egroup + \par + } +%% <<< history <<< + +\let\cny\@firstofone +\let\dept\@firstofone +\let\institution\@firstofone +\let\addr\@firstofone +\let\city\@firstofone + +\def\dedicated#1{% + \gdef\@dedicated{#1}% + } +\def\dedicated@fmt#1{% + \normalsize\itshape + #1% + } +%% <<< Frontmatter <<< + +%% fix footnote Postscript /Dest position +\let\orig@H@@footnotetext\H@@footnotetext +\long\def\H@@footnotetext#1{{% + \def\@pdfview{XYZ}% + \def\@pdfviewparams{ pdf@hoff pdf@voff \strip@pt\baselineskip\space 72.27 div 72 mul add null}% + \orig@H@@footnotetext{#1}% + }} + +%% definitions execution after \begin{document} (for tex4ht) +\ifdefined\HCode + { + \catcode`#=12\relax + \catcode`&=6\relax + \gdef\startlocaldefs{% + \catcode`#=12\relax + \makeatletter + \@startlocaldefs + } + \gdef\@startlocaldefs@toks{}% + \long\gdef\@startlocaldefs&1\endlocaldefs{% + \long\xdef\@startlocaldefs@toks{\detokenize{&1}}% + \makeatother + \catcode`#=6\relax + } + } + \appto\document{\makeatletter\expandafter\scantokens\expandafter{\@startlocaldefs@toks}\makeatother}% +\fi + +%% Acknowledgement +\define@key{acknowledgement}{title}{\def\acknowledgement@title{#1}} +\def\acknowledgement@title{Acknowledgement} +\newenvironment{acknowledgement}[1][]% + {% + \setkeys{acknowledgement}{#1}% + \section*{\acknowledgement@title}% + }{% + \par + } + +%% >>> Funding >>> +\def\funding@font{\normalfont \Large \bfseries} + +\def\@afterbox{% + \everypar{% + \if@nobreak + \@nobreakfalse + \clubpenalty \@M + \if@afterindent \else + {\setbox\z@\lastbox}% + \everypar{}% + \fi + \else + \clubpenalty \@clubpenalty + {\setbox\z@\lastbox}% + \everypar{}% + \fi + }% + } + +\newenvironment{Env@EndOfContrib}[1] + {% + \if@noskipsec\leavevmode\fi + \par + \@afterindenttrue + \if@nobreak + \everypar{}% + \else + \addpenalty\@secpenalty\addvspace\baselineskip + \fi + \@nobreakfalse + \global\@noskipsectrue + \everypar{% + \if@noskipsec + \global\@noskipsecfalse + {\setbox\z@\lastbox}% + \clubpenalty\@M + {% + %\sffamily + \funding@font + #1\\\\% + }% + %\unskip\hskip0.5em\relax + \else + \clubpenalty \@clubpenalty + \everypar{}% + \fi + }% + \ignorespaces + }{% + \par + \vskip\baselineskip + \aftergroup\@afterbox + } + +\define@key{gsponsor}{id}{\def\gsponsor@id{#1}} +\define@key{gsponsor}{sponsor-id}{\def\gsponsor@sponsorid{#1}} +\define@key{gnumber}{refid}{\def\gnumber@refid{#1}} + +\def\gsponsor@fmt#1{#1} +\def\gnumber@fmt#1{#1} + +\newcommand\gsponsor[2][]{% + \bgroup + \setkeys{gsponsor}{#1}% + \gsponsor@fmt{#2}% + \egroup + } +\newcommand\gnumber[2][]{% + \bgroup + \setkeys{gnumber}{#1}% + \gnumber@fmt{#2}% + \egroup + } + +\def\fundingname{Funding} +\def\funding{\Env@EndOfContrib\fundingname} +\let\endfunding\endEnv@EndOfContrib +%% <<< Funding <<< + +%% >>> esm environment for xml >>> +\def\esm@titlesep{} +\def\esm{% + \@ifnextchar[% + \@esm + {\@esm[]}% + } +\def\@esm[#1]{% + \bgroup + \par + \def\esmtitle##1{% + \section*{##1}% + \esm@titlesep}% + \def\esmdescription{\@firstofone}% + \def\esmfiletype{\@gobble}% + \def\esmfilename{\@gobble}% + } +\def\endesm{\egroup} +%% <<< esm environment for xml <<< + +\newenvironment{backmatter}{\goodbreak\ignorespaces}{} + +\newenvironment{biography}{% + \par\smallskip + \def\author##1{\noindent\textbf{##1}}% + \def\photofilename##1{% + \includegraphics[scale=1]{##1}\\ + }% + }% + {}% + +\newdimen\tablewidth +\tablewidth=\textwidth +\belowcaptionskip=4pt + +%%%%%%%%%%%%%%% CONTTABLE +\newtoks\table@head + +\long\def\thead#1\endthead{\global\table@head{#1}#1} +\def\contthead{\the\table@head} + +\let\contpart\relax + +\newenvironment{conttable}{% + \let\footnoterule\relax + \tablewidth=\hsize + \@float{table}} + {\end@float} + +\newenvironment{conttable*}{% + \let\footnoterule\relax + \tablewidth=\textwidth + \@dblfloat{table}} + {\end@dblfloat} + +\newenvironment{sidewaysconttable}{% + \@rotfloat{table}} + {\end@rotfloat} + +%%%%%%%%%%%%%%% CONTCAPTION + +\def\contcaptionbraces#1#2{\gdef\contcaptionbraces@left{#1}\gdef\contcaptionbraces@right{#2}} +\contcaptionbraces() + +\def\contcaption#1{\caption{\contcaptionbraces@left\textit{#1}\contcaptionbraces@right}} + +\def\continued#1{\rightline{(\textit{#1})}} + +\@namedef{contfigure*}{\@nameuse{figure*}} +\@namedef{endcontfigure*}{\@nameuse{endfigure*}} +\let\contfigure\figure +\let\endcontfigure\endfigure + +\let\tsub\textsubscript +\let\tsup\textsuperscript + +\RequirePackage{marginnote} +% \reversemarginpar +\setlength{\marginparwidth}{2.4cm} + +\endinput +%% +%% End of file `jds.cls'. diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/shortcodes.lua b/blockedcholesky/_extensions/wenjie2wang/jds/shortcodes.lua new file mode 100644 index 0000000..9bcb925 --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/shortcodes.lua @@ -0,0 +1,37 @@ +--[[ + Shortcodes for the Journal of Data Science Quarto format. + Provides cross-format support for common JDS markup. +]]-- + +-- Shortcode for nicely formatted 'LaTeX' string +function latex() + if quarto.doc.isFormat("pdf") then + return pandoc.RawInline('tex', '{\\LaTeX}') + elseif quarto.doc.isFormat("html") then + return pandoc.Math('InlineMath', "\\LaTeX") + else + return pandoc.Span('LaTeX') + end +end + +-- Shortcode for programming language names: {{< proglang R >}} +function proglang(args) + local lang = pandoc.utils.stringify(args[1]) + if lang == nil or lang == "" then return pandoc.Str("") end + if quarto.doc.isFormat("pdf") then + return pandoc.RawInline('tex', '\\proglang{' .. lang .. '}') + else + return pandoc.Strong(pandoc.Str(lang)) + end +end + +-- Shortcode for software package names: {{< pkg stats >}} +function pkg(args) + local name = pandoc.utils.stringify(args[1]) + if name == nil or name == "" then return pandoc.Str("") end + if quarto.doc.isFormat("pdf") then + return pandoc.RawInline('tex', '\\pkg{' .. name .. '}') + else + return pandoc.Emph(pandoc.Str(name)) + end +end diff --git a/blockedcholesky/_extensions/wenjie2wang/jds/template.tex b/blockedcholesky/_extensions/wenjie2wang/jds/template.tex new file mode 100644 index 0000000..1c7b7c0 --- /dev/null +++ b/blockedcholesky/_extensions/wenjie2wang/jds/template.tex @@ -0,0 +1,203 @@ +%% Custom Pandoc template for Journal of Data Science (jdsart.cls) +%% This template is used instead of Quarto's default because jdsart.cls +%% has its own frontmatter structure for title, authors, and affiliations. + +\documentclass[$for(classoption)$$classoption$$sep$,$endfor$]{jdsart} + +%% ----- JDS article metadata ----- +$if(firstpage)$ +\setcounter{page}{$firstpage$} +$else$ +\setcounter{page}{1} +$endif$ +$if(month)$ +\pubmonth{$month$} +$else$ +\pubmonth{} +$endif$ +$if(year)$ +\pubyear{$year$} +$else$ +\pubyear{20xx} +$endif$ +$if(volume)$ +\volume{$volume$} +$else$ +\volume{xx} +$endif$ +$if(issue)$ +\issue{$issue$} +$else$ +\issue{xx} +$endif$ +$if(doi)$ +\doi{$doi$} +$else$ +\doi{10.6339/xx-JDSxxxx} +$endif$ +$if(received)$\received{$received$}$endif$ +$if(revised)$\revised{$revised$}$endif$ +$if(accepted)$\accepted{$accepted$}$endif$ +$if(publishedonline)$\publishedonline{$publishedonline$}$endif$ + +%% ----- Syntax highlighting ----- +$if(highlighting-macros)$ +$highlighting-macros$ +$endif$ + +%% ----- Core packages ----- +%% Note: jdsart.cls already loads fontenc, textcomp, graphicx, xcolor, +%% natbib, hyperref, geometry, and fancyhdr. +$if(pdftex)$ +\usepackage[utf8]{inputenc} +$endif$ + +\providecommand{\tightlist}{% + \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}} + +$if(strikeout)$ +\usepackage[normalem]{ulem} +$endif$ + +%% ----- Graphics scaling (Quarto convention) ----- +$if(graphics)$ +\makeatletter +\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi} +\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi} +\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio} +\newsavebox\pandoc@box +\newcommand*\pandocbounded[1]{% scales image to fit in text height/width + \sbox\pandoc@box{#1}% + \Gscale@div\@tempa{\textheight}{\dimexpr\ht\pandoc@box+\dp\pandoc@box\relax}% + \Gscale@div\@tempb{\linewidth}{\wd\pandoc@box}% + \ifdim\@tempb\p@<\@tempa\p@\let\@tempa\@tempb\fi% + \ifdim\@tempa\p@<\p@\scalebox{\@tempa}{\usebox\pandoc@box}% + \else\usebox{\pandoc@box}% + \fi% +} +\makeatother +$endif$ + +%% ----- Table support ----- +$if(tables)$ +\usepackage{longtable,booktabs,array} +$if(multirow)$ +\usepackage{multirow} +$endif$ +\usepackage{calc} +\usepackage{footnote} +\makesavenoteenv{longtable} +$endif$ + +%% ----- CSL references (non-natbib fallback) ----- +$if(csl-refs)$ +\newlength{\cslhangindent} +\setlength{\cslhangindent}{1.5em} +\newlength{\csllabelwidth} +\setlength{\csllabelwidth}{3em} +\newenvironment{CSLReferences}[2] + {\begin{list}{}% + {\setlength{\itemindent}{0pt} + \setlength{\leftmargin}{0pt} + \ifnum#1>0 + \setlength{\leftmargin}{\cslhangindent} + \setlength{\itemindent}{-1\cslhangindent} + \fi + \setlength{\parsep}{\parskip} + \setlength{\itemsep}{#2\baselineskip}}} + {\end{list}} +$endif$ + +%% ----- Header includes from YAML ----- +$for(header-includes)$ +$header-includes$ +$endfor$ + +%% Reset Quarto's default ruled float style for code listings +\makeatletter +\@ifundefined{c@codelisting}{}{\floatstyle{plain}\restylefloat{codelisting}} +\makeatother + +%% ----- User preamble from YAML ----- +$if(preamble)$ +$preamble$ +$endif$ + +%% ===== Document begins ===== +\begin{document} +\begin{frontmatter} + +$if(pretitle)$ +\pretitle{$pretitle$} +$endif$ + +\title{$title$} +$if(short-title)$\runtitle{$short-title$}$endif$ + +$for(author)$ +\author[$author.marker$]{% + \inits{$author.first-init$}% + \fnms{$author.first$}% + \snm{$author.last$}% +$if(author.footnote-id)$ + \thanksref{$author.footnote-id$}% +$endif$ +$if(author.email)$ + \ead{$author.email$}% +$endif$ +} +$endfor$ + +$for(footnote)$ +\thankstext[$if(footnote.type)$type=$footnote.type$,$endif$id=$footnote.id$]{$footnote.content$} +$endfor$ + +$for(affiliation)$ +\address[$affiliation.marker$]{% +$if(affiliation.prefix)$$affiliation.prefix$, $endif$% + \institution{$affiliation.institution$}, \cny{$affiliation.country$}} +$endfor$ + +$if(abstract)$ +\begin{abstract} +$abstract$ +\end{abstract} +$endif$ + +$if(keywords)$ +\begin{keywords} +$for(keywords)$\kwd{$keywords$}$endfor$. +\end{keywords} +$endif$ + +\end{frontmatter} + +$for(include-before)$ +$include-before$ +$endfor$ + +$body$ + +$if(natbib)$ +$if(bibliography)$ +$if(biblio-title)$ +\renewcommand\refname{$biblio-title$} +$endif$ +$if(biblio-style)$ +\bibliographystyle{$biblio-style$} +$else$ +\bibliographystyle{jds} +$endif$ +\bibliography{$for(bibliography)$$bibliography$$sep$,$endfor$} +$endif$ +$endif$ + +$if(biblatex)$ +\printbibliography$if(biblio-title)$[title=$biblio-title$]$endif$ +$endif$ + +$for(include-after)$ +$include-after$ +$endfor$ + +\end{document} diff --git a/blockedcholesky/bibliography.bib b/blockedcholesky/bibliography.bib index c131f39..1f4b0e3 100644 --- a/blockedcholesky/bibliography.bib +++ b/blockedcholesky/bibliography.bib @@ -555,7 +555,7 @@ @article{JSSv107i04 @article{10.1145/1391989.1391995, author = {Chen, Yanqing and Davis, Timothy A. and Hager, William W. and Rajamanickam, Sivasankaran}, -title = {Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate}, +title = {Algorithm 887: {CHOLMOD}, Supernodal Sparse {C}holesky Factorization and Update/Downdate}, year = {2008}, issue_date = {October 2008}, publisher = {Association for Computing Machinery}, diff --git a/blockedcholesky/data/config.txt b/blockedcholesky/data/config.txt new file mode 100644 index 0000000..2b7daee --- /dev/null +++ b/blockedcholesky/data/config.txt @@ -0,0 +1,48 @@ +julia> BLAS.get_config() +LinearAlgebra.BLAS.LBTConfig +Libraries: +├ [ILP64] libmkl_rt.so +└ [ LP64] libmkl_rt.so + +julia> versioninfo() +Julia Version 1.12.6 +Commit 15346901f00 (2026-04-09 19:20 UTC) +Build Info: + Official https://julialang.org release +Platform Info: + OS: Linux (x86_64-linux-gnu) + CPU: 32 × AMD EPYC-Milan Processor + WORD_SIZE: 64 + LLVM: libLLVM-18.1.7 (ORCJIT, znver3) + GC: Built with stock GC +Threads: 12 default, 1 interactive, 12 GC (on 32 virtual cores) + +julia> versioninfo(verbose=true) +Julia Version 1.12.6 +Commit 15346901f00 (2026-04-09 19:20 UTC) +Build Info: + Official https://julialang.org release +Platform Info: + OS: Linux (x86_64-linux-gnu) + Ubuntu 26.04 LTS + uname: Linux 7.0.0-15-generic #15-Ubuntu SMP PREEMPT_DYNAMIC Wed Apr 22 16:06:43 UTC 2026 x86_64 unknown + CPU: AMD EPYC-Milan Processor: + speed user nice sys idle irq + #1-32 0 MHz 1992427 s 11 s 3398 s 7530497 s 0 s + Memory: 122.82613754272461 GB (123986.328125 MB free) + Uptime: 297761.02 sec + Load Avg: 0.02 0.01 0.0 + WORD_SIZE: 64 + LLVM: libLLVM-18.1.7 (ORCJIT, znver3) + GC: Built with stock GC +Threads: 12 default, 1 interactive, 12 GC (on 32 virtual cores) + +(ml32m) pkg> st +Status `~/ml32m/Project.toml` + [69666777] Arrow v2.8.1 + [336ed68f] CSV v0.10.16 + [a93c6f00] DataFrames v1.8.2 + [ff71e718] MixedModels v5.5.1 + [7e9fb7ac] MixedModelsDatasets v0.1.2 + [a5390f91] ZipFile v0.10.1 + [f43a241f] Downloads v1.7.0 diff --git a/blockedcholesky/data/moviesnratngs.arrow b/blockedcholesky/data/moviesnratngs.arrow index 658f2b7..e9d799d 100644 Binary files a/blockedcholesky/data/moviesnratngs.arrow and b/blockedcholesky/data/moviesnratngs.arrow differ diff --git a/blockedcholesky/data/sizespeed.csv b/blockedcholesky/data/sizespeed.csv index 1654151..40aeb8e 100644 --- a/blockedcholesky/data/sizespeed.csv +++ b/blockedcholesky/data/sizespeed.csv @@ -1,22 +1,22 @@ -mc,uc,nratings,nusers,nmvie,modelsz,L22sz,nv,fittime,evtime -1,20,32000204,200948,84432,56.24513,53.113422,26,18740.139,667.45874 -2,20,31981597,200948,65825,35.411522,32.28285,19,7091.5474,355.21213 -5,20,31921467,200948,43884,17.469769,14.348369,26,3816.2913,142.36954 -10,20,31842705,200948,31961,10.7227545,7.6108093,25,2165.8496,82.710724 -15,20,31777786,200948,26428,8.30907,5.2037773,24,1619.6187,64.58687 -20,20,31725920,200948,23350,7.1622887,4.0622244,21,1253.507,57.18353 -50,20,31498689,200947,16034,4.992517,1.9154636,28,1212.5883,42.08827 -1,40,30433400,144848,84205,55.799446,52.82821,23,15952.788,663.37054 -2,40,30415014,144848,65819,35.24521,32.276962,32,11565.444,356.29413 -5,40,30355623,144848,43884,17.30941,14.348369,38,5496.8467,139.71104 -10,40,30277758,144848,31961,10.562482,7.6108093,21,1840.2719,83.437744 -15,40,30213583,144848,26428,8.148869,5.2037773,21,1428.1853,64.851364 -20,40,30162330,144848,23350,7.0021467,4.0622244,17,1026.0945,57.07803 -50,40,29938038,144848,16034,4.8326597,1.9154636,24,1046.2183,41.569332 -1,80,27569316,94380,83897,55.13217,52.44245,23,16155.171,677.0729 -2,80,27551218,94380,65799,34.944122,32.25735,22,8541.434,392.76468 -5,80,27492886,94380,43884,17.028042,14.348369,23,3564.9575,142.58371 -10,80,27416310,94380,31961,10.2812395,7.6108093,23,2054.1672,85.285385 -15,80,27353366,94380,26428,7.867746,5.2037773,26,1785.077,65.71991 -20,80,27303019,94380,23350,6.721114,4.0622244,43,2527.2083,57.717804 -50,80,27083536,94380,16034,4.55209,1.9154636,20,873.97784,41.380127 +mc,uc,nratings,nusers,nmvie,modelsz,L22sz,nv,fittime,mnevtm +1,20,32000204,200948,84432,55.999664,53.113422,30,9129.769,304.32562 +2,20,31981597,200948,65825,35.1666,32.28285,39,6546.0786,167.84818 +5,20,31921467,200948,43884,17.225552,14.348369,32,2269.9685,70.936516 +10,20,31842705,200948,31961,10.479625,7.6108093,33,1484.9624,44.99886 +15,20,31777786,200948,26428,8.066437,5.2037773,35,1248.0991,35.659973 +20,20,31725920,200948,23350,6.920109,4.0622244,34,1106.2368,32.536377 +50,20,31498689,200947,16034,4.7523084,1.9154636,35,885.96277,25.313223 +1,40,30433400,144848,84205,55.567394,52.82821,35,10657.067,304.48764 +2,40,30415014,144848,65819,35.013683,32.276962,33,5474.281,165.8873 +5,40,30355623,144848,43884,17.07859,14.348369,39,2789.3848,71.52269 +10,40,30277758,144848,31961,10.332743,7.6108093,37,1707.9968,46.162075 +15,40,30213583,144848,26428,7.919621,5.2037773,29,1049.0757,36.175022 +20,40,30162330,144848,23350,6.773348,4.0622244,32,1048.7473,32.773354 +50,40,29938038,144848,16034,4.60581,1.9154636,37,900.5709,24.339754 +1,80,27569316,94380,83897,54.92231,52.44245,30,9101.557,303.38522 +2,80,27551218,94380,65799,34.734783,32.25735,31,5126.242,165.36266 +5,80,27492886,94380,43884,16.819399,14.348369,32,2209.0645,69.033264 +10,80,27416310,94380,31961,10.073666,7.6108093,39,1729.2856,44.340656 +15,80,27353366,94380,26428,7.6606545,5.2037773,30,1100.8684,36.695614 +20,80,27303019,94380,23350,6.5144625,4.0622244,33,1076.3864,32.617767 +50,80,27083536,94380,16034,4.3473544,1.9154636,30,717.3222,23.91074 diff --git a/blockedcholesky/data/usersnratngs.arrow b/blockedcholesky/data/usersnratngs.arrow index d6b38ca..afd8656 100644 Binary files a/blockedcholesky/data/usersnratngs.arrow and b/blockedcholesky/data/usersnratngs.arrow differ diff --git a/blockedcholesky/replication/benchmark.R b/blockedcholesky/replication/benchmark.R deleted file mode 100644 index 9a8a8c7..0000000 --- a/blockedcholesky/replication/benchmark.R +++ /dev/null @@ -1,19 +0,0 @@ -library(lme4, quietly = TRUE) -library(glmmTMB) -library(microbenchmark) -dat <- InstEval -dat$service <- as.numeric(dat$service) - 1 # convert to a 0/1 numeric vector -ctrl <- lmerControl(calc.derivs = FALSE) -form <- y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) -mlmer <- lmer(form, dat, REML = FALSE, control = ctrl) -print(summary(mlmer, correlation = FALSE)) -tmbmod <- glmmTMB(form, dat) -print(summary(tmbmod)) -microbenchmark( - lmer = lmer(form, dat, REML = FALSE, control = ctrl), - glmmTMB = glmmTMB(form, dat), times = 6 -) |> print() -# R and package versions -print(R.version.string) -print(packageVersion("lme4")) -print(packageVersion("glmmTMB")) \ No newline at end of file diff --git a/blockedcholesky/replication/replication.jl b/blockedcholesky/replication/replication.jl deleted file mode 100644 index b5adebe..0000000 --- a/blockedcholesky/replication/replication.jl +++ /dev/null @@ -1,379 +0,0 @@ -using Pkg -Pkg.instantiate() - -Pkg.add("RCall") - - -#| output: false -using AlgebraOfGraphics, CairoMakie, PrettyTables, Printf, SparseArrays, LinearAlgebra -using Arrow, Chairmarks, CSV, DataFrames, MixedModels -dat = DataFrame(MixedModels.dataset("insteval")) -dat.service = float.(dat.service .== "Y") # convert to 0/1 encoding -sumry = describe(dat, :min, :max, :mean, :nunique, :eltype) - -#| echo: false -#| output: asis -#| tbl-cap: Summary of the `insteval` dataset -#| label: tbl-insteval-summary -fmt_nothing = (v, i, j) -> isnothing(v) ? "-" : v -pretty_table( - sumry, - formatters = (ft_round(2, [4]), fmt_nothing), - show_subheader = false, backend = Val(:latex)) - -#| echo: false -#| output: false -# set up the environment -using MKL_jll -hostarch = Base.BinaryPlatforms.arch(Base.BinaryPlatforms.HostPlatform()) -@static if Sys.isapple() && hostarch == "aarch64" - using AppleAccelerate -elseif MKL_jll.is_available() - using MKL -end -# override QuartoNotebookRunner https://github.com/PumasAI/QuartoNotebookRunner.jl/issues/294 -CairoMakie.activate!(; type="pdf") - -form = @formula(y ~ 1 + service + (1|d) + (1|s) + (1|dept) + (0 + service|dept)) -m1 = fit(MixedModel, form, dat, progress=false) # suppress display of a progress bar -print(m1) - -@be fit($MixedModel, $form, $dat; progress=false) seconds=8 - -@be objective(updateL!(setθ!($m1, $(m1.θ)))) - -show([typeof(rt) for rt in m1.reterms]) - -show([length(rt.levels) for rt in m1.reterms]) - -show([size(rt, 2) for rt in m1.reterms]) - -#| tbl-cap: Block Structure of the `A` matrix \label{tab:blk-descrip} -#| tbl-cap-location: top -BlockDescription(m1) - -θ̂ = m1.θ -show(θ̂) - -m1.λ - -#| echo: false -#| label: fit-false -m2 = LinearMixedModel(form, dat); - -#| echo: false -#| label: fit-reorder -m2.reterms[1:2] = m2.reterms[2:-1:1]; # swap elements 1 and 2 - -A, L = MixedModels.createAL(m2.reterms, m2.Xymat); # recreate A, L - -copyto!(m2.A, A); # set A -copyto!(m2.L, L); # set L - -#| echo: false -refit!(m2; progress = false); - -#| echo: false -#| fig-cap: Sparsity pattern of the augmented matrix $\mbfA$ when the order of the random effects is not carefully chosen (left) and when the order is carefully chosen (right). -#| label: fig-a -include("scripts/makiespy.jl") -CairoMakie.activate!(; type="png") - -f = Figure() -spyA(m2, f[1,1]; markersize = 5, colormap = :imola) # m2 is the slower model -spyA(m1, f[1,2]; markersize = 5, colormap = :imola) # m1 is the slower model -f - -#| echo: false -#| fig-cap: Sparsity pattern of the lower triangular Cholesky factor when the order of the random effects is not carefully chosen (left) and when the order is carefully chosen (right). -#| label: fig-blk - -second(x) = x[2] -reszfast = second.(size.(m1.reterms)) -reszslow = second.(size.(m2.reterms)) -szl, _ = size(sparseL(m1; full = true)) - - -f = Figure(figure_padding = 60) -spyL(m2, f[1,1]; markersize = 5, colormap = :imola) # m2 is the slower model -spyL(m1, f[1,2]; markersize = 5, colormap = :imola) # m1 is the slower model -colgap!(f.layout, 50.) - - -x0, y0 = content(f[1,1]).scene.viewport.val.origin -xw, yw = content(f[1,1]).scene.viewport.val.widths -bracket!( - f.scene, - x0, y0 + yw, - x0 + xw * (reszslow[1] / szl), y0 + yw, - text = string(reszslow[1]) -) - -bracket!( - f.scene, - x0 + xw * (reszslow[1] / szl), y0 + yw, - x0 + xw * (reszslow[1] + reszslow[2]) / szl, y0 + yw, - text = string(reszslow[2]) -) - - - -bracket!( - f.scene, - x0, y0 + yw, - x0, y0 + yw - yw * (reszslow[1] / szl), - text = string(reszslow[1]), - orientation = :down, -) - -bracket!( - f.scene, - x0, y0 + yw - yw * (reszslow[1] / szl), - x0, y0 + yw - yw * (reszslow[1] + reszslow[2]) / szl, - text = string(reszslow[2]), - orientation = :down -) - -x0, y0 = content(f[1,2]).scene.viewport.val.origin -xw, yw = content(f[1,2]).scene.viewport.val.widths -bracket!( - f.scene, - x0, y0 + yw, - x0 + xw * reszfast[1] / szl, y0 + yw, - text = string(reszfast[1]) -) - -bracket!( - f.scene, - x0 + xw * reszfast[1] / szl, y0 + yw, - x0 + xw * (reszfast[1] + reszfast[2]) / szl, y0 + yw, - text = string(reszfast[2]) -) - -bracket!( - f.scene, - x0, y0 + yw, - x0, y0 + yw - yw * reszfast[1] / szl , - text = string(reszfast[1]), - orientation = :down -) - -bracket!( - f.scene, - x0, y0 + yw - yw * reszfast[1] / szl, - x0, y0 + yw - yw * (reszfast[1] + reszfast[2]) / szl, - text = string(reszfast[2]), - orientation = :down -) - -f - -#| echo: false -CairoMakie.activate!(; type="pdf"); - -#| echo: false -#| output: false -movies = DataFrame(Arrow.Table("./data/moviesnratngs.arrow")) -users = DataFrame(Arrow.Table("./data/usersnratngs.arrow")) - -#| echo: false -#| fig-cap: "Empirical distribution plots of the number of ratings per movie and per user. The horizontal axes are on a logarithmic scale." -#| label: fig-nrtngsecdf -#| warning: false -let - f = Figure(; size=(800, 300)) - xscale = log10 - xminorticksvisible = true - xminorgridvisible = true - yminorticksvisible = true - xminorticks = IntervalsBetween(10) - ylabel = "Relative cumulative frequency" - nrtngs = sort(movies.nrtngs) - ecdfplot( - f[1, 1], - nrtngs; - npoints=last(nrtngs), - axis=( - xlabel="Number of ratings per movie (logarithmic scale)", - xminorgridvisible, - xminorticks, - xminorticksvisible, - xscale, - ylabel, - yminorticksvisible, - ), - ) - urtngs = sort(users.nrtngs) - ecdfplot( - f[1, 2], - urtngs; - npoints=last(urtngs), - axis=( - xlabel="Number of ratings per user (logarithmic scale)", - xminorgridvisible, - xminorticks, - xminorticksvisible, - xscale, - yminorticksvisible, - ), - ) - f -end - -#| echo: false -#| output: false -sizespeed = CSV.read("./data/sizespeed.csv", DataFrame; downcast=true); - -#| echo: false -#| label: fig-nratingsbycutoff -#| fig-cap: "Number of ratings (left) and movies (right) in reduced table by movie cutoff and by user cutoff" -#| warning: false - -function siformat(x) - if x >= 1_000_000 - val = x / 1_000_000 - suffix = "M" - elseif x >= 1_000 - val = x / 1_000 - suffix = "K" - else - val = x - suffix = "" - end - - return string(@sprintf("%g", val), suffix) -end - -draw( - data(sizespeed) * - mapping( - :mc => "Minimum number of ratings per movie", - [:nratings => "Total number of ratings", :nmvie => "Number of movies in table"], - col = dims(1) => renamer(["Ratings", "Movies"]), - color = :uc => renamer([20 => "≥20 ratings", 40 => "≥40 ratings", 80 => "≥80 ratings"]) => "per user", - ) * visual(ScatterLines); - legend = (; position = :bottom, titleposition = :left), - axis = (; ytickformat = values -> [siformat(value) for value in values]), - figure=(; size=(800, 350)) -) - -#| echo: false -#| fig-cap: Memory footprint of the model representation by minimum number of ratings per user and per movie. -#| label: fig-memoryfootprint -#| warning: false -draw( - data(sizespeed) * - mapping( - :mc => "Minimum number of ratings per movie", - :modelsz => "Size of model object (GiB)", - color = :uc => renamer([20 => "≥20 ratings", 40 => "≥40 ratings", 80 => "≥80 ratings"]) => "per user",) * - visual(ScatterLines); - legend = (; position = :bottom, titleposition = :left), -) - -#| echo: false -#| fig-cap: Proportion of memory footprint of the model in L[2,2] versus the overall model size (GiB). -#| label: fig-l22prop -#| warning: false -draw( - data(transform!(sizespeed, [:L22sz, :modelsz] => ((x, y) -> x ./ y) => :L22prop)) * - mapping( - :modelsz => "Size of model object (GiB)", - :L22prop => "Proportion of memory footprint in L[2,2]", - color = :uc => renamer([20 => "≥20 ratings", 40 => "≥40 ratings", 80 => "≥80 ratings"]) => "per user",) * - visual(ScatterLines); - legend = (; position = :bottom, titleposition = :left), -) - -#| echo: false -#| fig-cap: "Evaluation time for the objective (s) versus size of the [2,2] block of L (GiB)" -#| label: fig-evtimevsl22 -#| warning: false -draw( - data(sizespeed) * - mapping( - :L22sz => "Size of [2,2] block of L (GiB)", - :evtime => "Time for one evaluation of objective (s)", - color = :uc => renamer([20 => "≥20 ratings", 40 => "≥40 ratings", 80 => "≥80 ratings"]) => "per user",) * - visual(ScatterLines); - legend = (; position = :bottom, titleposition = :left), -) - -mreml = fit( - MixedModel, - @formula(y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)), - dat, - REML = true, - progress=false -) -print(mreml) - -#| eval: false -m2 = fit( - MixedModel, - @formula(y ~ 1 + service + (1 | d) + (1 | s) + (1 | dept) + (0 + service | dept)), - dat, - progress=false # suppress the display of a progress bar -); - -#| eval: false -m2.reterms[1:2] = m2.reterms[2:-1:1]; # swap elements 1 and 2 - -A, L = MixedModels.createAL(m2.reterms, m2.Xymat); # recreate A, L - -copyto!(m2.A, A); # set A -copyto!(m2.L, L); # set L - -@be refit!($m2; progress = false) seconds = 100 # very slow - -@be refit!($m1, progress = false) seconds = 8 - -print(m2) - -#| echo: false -#| output: asis -#| tbl-cap: Summary of results for the `ml-32m` data modeling run speeds -#| label: tbl-sizespeed -sizespeed = CSV.read("./data/sizespeed.csv", DataFrame; downcast=true); -rename!(sizespeed, - :mc => "movie cutoff", - :uc => "user cutoff", - :nratings => "ratings", - :nusers => "users", - :nmvie => "movies", - :modelsz => "model (GiB)", - :L22sz => "L[2,2] (GiB)", - :fittime => "time (s)", - :nv => "n eval", - :evtime => "time per eval (s)") -pretty_table(sizespeed; formatters = ft_round(2, [6, 7, 9, 10]), backend = Val(:markdown), show_subheader = false) - - - -BLAS.get_config() - -using Pkg; Pkg.status() - -using RCall - -R""" -library(lme4, quietly = TRUE) -library(glmmTMB) -library(microbenchmark) -dat <- InstEval -dat$service <- as.numeric(dat$service) - 1 # convert to a 0/1 numeric vector -ctrl <- lmerControl(calc.derivs = FALSE) -form <- y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) -mlmer <- lmer(form, dat, REML = FALSE, control = ctrl) -print(summary(mlmer, correlation = FALSE)) -tmbmod <- glmmTMB(form, dat) -print(summary(tmbmod)) -microbenchmark( - lmer = lmer(form, dat, REML = FALSE, control = ctrl), - glmmTMB = glmmTMB(form, dat), times = 6 -) |> print() -# R and package versions -print(R.version.string) -print(packageVersion("lme4")) -print(packageVersion("glmmTMB")) -""" diff --git a/blockedcholesky/scripts/benchmark.R b/blockedcholesky/scripts/benchmark.R deleted file mode 100644 index 11c058d..0000000 --- a/blockedcholesky/scripts/benchmark.R +++ /dev/null @@ -1,19 +0,0 @@ -library(lme4, quietly = TRUE) -library(glmmTMB) -library(microbenchmark) -dat <- InstEval -dat$service <- as.numeric(dat$service) - 1 # convert to a 0/1 numeric vector -ctrl <- lmerControl(calc.derivs = FALSE) -form <- y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) -mlmer <- lmer(form, dat, REML = FALSE, control = ctrl) -summary(mlmer, correlation = FALSE) -tmbmod <- glmmTMB(form, dat) -summary(tmbmod) -microbenchmark( - lmer = lmer(form, dat, REML = FALSE, control = ctrl), - glmmTMB = glmmTMB(form, dat), times = 6 -) -# R and package versions -R.version.string -packageVersion("lme4") -packageVersion("glmmTMB") \ No newline at end of file diff --git a/blockedcholesky/scripts/benchmarkoutput.txt b/blockedcholesky/scripts/benchmarkoutput.txt deleted file mode 100644 index 65de969..0000000 --- a/blockedcholesky/scripts/benchmarkoutput.txt +++ /dev/null @@ -1,80 +0,0 @@ -> library(lme4, quietly = TRUE) -> library(glmmTMB) -> library(microbenchmark) -> dat <- InstEval -> dat$service <- as.numeric(dat$service) - 1 # convert to a 0/1 numeric vector -> ctrl <- lmerControl(calc.derivs = FALSE) -> form <- y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) -> mlmer <- lmer(form, dat, REML = FALSE, control = ctrl) -> summary(mlmer, correlation = FALSE) -Linear mixed model fit by maximum likelihood ['lmerMod'] -Formula: y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | - dept) - Data: dat -Control: ctrl - - AIC BIC logLik -2*log(L) df.resid - 237662.6 237727.0 -118824.3 237648.6 73414 - -Scaled residuals: - Min 1Q Median 3Q Max --2.99959 -0.74769 0.04009 0.77283 3.11345 - -Random effects: - Groups Name Variance Std.Dev. - s (Intercept) 0.105295 0.32449 - d (Intercept) 0.262425 0.51227 - dept (Intercept) 0.002579 0.05079 - dept.1 service 0.023405 0.15299 - Residual 1.385009 1.17686 -Number of obs: 73421, groups: s, 2972; d, 1128; dept, 14 - -Fixed effects: - Estimate Std. Error t value -(Intercept) 3.27765 0.02350 139.463 -service -0.05074 0.04399 -1.153 -> tmbmod <- glmmTMB(form, dat) -> summary(tmbmod) - Family: gaussian ( identity ) -Formula: -y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept) + (0 + service | dept) -Data: dat - - AIC BIC logLik -2*log(L) df.resid - 237662.6 237727.0 -118824.3 237648.6 73414 - -Random effects: - -Conditional model: - Groups Name Variance Std.Dev. - s (Intercept) 0.10530 0.32449 - d (Intercept) 0.26243 0.51228 - dept (Intercept) 0.00258 0.05079 - dept.1 service 0.02340 0.15296 - Residual 1.38501 1.17686 -Number of obs: 73421, groups: s, 2972; d, 1128; dept, 14 - -Dispersion estimate for gaussian family (sigma^2): 1.39 - -Conditional model: - Estimate Std. Error z value Pr(>|z|) -(Intercept) 3.27765 0.02351 139.39 <2e-16 *** -service -0.05074 0.04408 -1.15 0.25 ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -> microbenchmark( -+ lmer = lmer(form, dat, REML = FALSE, control = ctrl), -+ glmmTMB = glmmTMB(form, dat), times = 6 -+ ) -Unit: seconds - expr min lq mean median uq max neval cld - lmer 17.61183 17.61670 17.81537 17.78631 18.01708 18.07398 6 a - glmmTMB 23.38106 23.47987 23.56080 23.52521 23.62934 23.82409 6 b -> # R and package versions -> R.version.string -[1] "R version 4.5.0 (2025-04-11)" -> packageVersion("lme4") -[1] ‘1.1.37’ -> packageVersion("glmmTMB") -[1] ‘1.1.11’ -> diff --git a/blockedcholesky/scripts/ml32-download.jl b/blockedcholesky/scripts/ml32-download.jl new file mode 100644 index 0000000..2556dd3 --- /dev/null +++ b/blockedcholesky/scripts/ml32-download.jl @@ -0,0 +1,34 @@ +using Arrow, CSV, DataFrames, Downloads, MixedModels, ZipFile +const ML_32M_URL = "https://files.grouplens.org/datasets/movielens/ml-32m.zip" + +const ML_32M_LOCAL = joinpath("data", "ml-32m.zip") + +function extract_csv(zipfile, fname; kwargs...) + file = only(filter(f -> endswith(f.name, fname), zipfile.files)) + return CSV.read(file, DataFrame; delim=',', header=1, kwargs...) +end + +function load_zipfile() + if !isfile(ML_32M_LOCAL) + @info "Downloading data" + Downloads.download(ML_32M_URL, ML_32M_LOCAL) + end + zipfile = ZipFile.Reader(ML_32M_LOCAL) + ratings = extract_csv( + zipfile, + "ratings.csv"; + drop=[4], + types=[Int32, Int32, Float32, Int32], + pool=[true, true, true, false], + ) + close(zipfile) + movies = combine(groupby(ratings, :movieId), nrow => :mnrtngs) # number of ratings per movie + movies.mnrtngs = Int32.(movies.mnrtngs) + users = combine(groupby(ratings, :userId), nrow => :unrtngs) + users.unrtngs = Int32.(users.unrtngs) + disallowmissing!(leftjoin!(leftjoin!(ratings, movies; on=:movieId), users; on=:userId)) + Arrow.write("data/ratings.arrow", ratings; compress=:zstd) + Arrow.write("data/moviesnratngs.arrow", movies; compress=:zstd) + Arrow.write("data/usersnratngs.arrow", users; compress=:zstd) + return ratings, movies, users +end diff --git a/blockedcholesky/scripts/ml32-timefit.jl b/blockedcholesky/scripts/ml32-timefit.jl new file mode 100644 index 0000000..aee0e71 --- /dev/null +++ b/blockedcholesky/scripts/ml32-timefit.jl @@ -0,0 +1,60 @@ +using Arrow, CSV, DataFrames, MixedModels, MKL + +function timefit( + mc::Integer, + uc::Integer, + ratings::AbstractDataFrame; + initial_step::Vector{<:AbstractFloat} = [0.5, 0.5], + ftol_rel::AbstractFloat = 1.0e-10, + backend::Symbol = :nlopt, + optimizer::Symbol = :LN_BOBYQA, + ) + mc = Int8(mc) + uc = Int8(uc) + df = ratings[(ratings.mnrtngs .≥ mc) .& (ratings.unrtngs .≥ uc), :] + nratings = Int32(size(df, 1)) + nusers = Int32(length(unique(df.userId))) + nmvie = Int32(length(unique(df.movieId))) + model = LinearMixedModel(@formula(rating ~ 1 + (1|userId) + (1|movieId)), df) + model.optsum.initial_step = initial_step + model.optsum.ftol_rel = ftol_rel + model.optsum.ftol_abs = 1.0e-6 + model.optsum.backend = backend + model.optsum.optimizer = optimizer + modelsz = Float32(Base.summarysize(model) / (2^30)) # size in GiB + L22sz = Float32(Base.summarysize(model.L[3]) / (2^30)) + fittime = Float32(@elapsed fit!(model; progress=isinteractive())) + nv = Int8(length(model.optsum.fitlog)) + mnevtm = fittime / nv + return (; mc, uc, nratings, nusers, nmvie, modelsz, L22sz, nv, fittime, mnevtm) +end + +const ratings = DataFrame(Arrow.Table("./data/ratings.arrow")) + +function mktbl() + res = @NamedTuple{ + mc::Int8, + uc::Int8, + nratings::Int32, + nusers::Int32, + nmvie::Int32, + modelsz::Float32, + L22sz::Float32, + nv::Int8, + fittime::Float32, + mnevtm::Float32, + }[] + for uc in Int8.([20, 40, 80]) + for mc in Int8.([1, 2, 5, 10, 15, 20, 50]) + try + sizespeedrow = timefit(mc, uc, ratings) + @info sizespeedrow + push!(res, sizespeedrow) + catch e + finally + CSV.write("./data/sizespeed.csv", res) + end + end + end +end + diff --git a/blockedcholesky/texincludes.tex b/blockedcholesky/texincludes.tex index 02ba325..aa20c80 100644 --- a/blockedcholesky/texincludes.tex +++ b/blockedcholesky/texincludes.tex @@ -41,11 +41,32 @@ \DeclareUnicodeCharacter{2303}{\^} \else \usepackage{fontspec} + \usepackage{bm} \setmonofont{JuliaMono} - \newcommand\mcN{{\mathbfcal{N}}} - \newcommand\mcB{{\mathbfcal{B}}} - \newcommand\mcU{{\mathbfcal{U}}} - \newcommand\mcY{{\mathbfcal{Y}}} + \newcommand\mbfZ{{\mathbf{Z}}} + \newcommand\mbfA{{\mathbf{A}}} + \newcommand\mbfI{{\mathbf{I}}} + \newcommand\mbfL{{\mathbf{L}}} + \newcommand\mbfO{{\mathbf{O}}} + \newcommand\mbfR{{\mathbf{R}}} + \newcommand\mbfS{{\mathbf{S}}} + \newcommand\mbfT{{\mathbf{T}}} + \newcommand\mbfX{{\mathbf{X}}} + \newcommand\mbfb{{\mathbf{b}}} + \newcommand\mbft{{\mathbf{t}}} + \newcommand\mbfu{{\mathbf{u}}} + \newcommand\mbfv{{\mathbf{v}}} + \newcommand\mbfy{{\mathbf{y}}} + \newcommand\mcN{\bm{\mathcal{N}}} + \newcommand\mcB{\bm{\mathcal{B}}} + \newcommand\mcU{\bm{\mathcal{U}}} + \newcommand\mcY{\bm{\mathcal{Y}}} + \newcommand\mbfLambda{{\bm{\Lambda}}} + \newcommand\mbfOmega{{\bm{\Omega}}} + \newcommand\mbfSigma{{\bm{\Sigma}}} + \newcommand\mbfbeta{{\bm{\beta}}} + \newcommand\mbftheta{{\bm{\theta}}} + \usepackage{newunicodechar} % linebreak added here to correct for https://github.com/JuliaLang/LinearAlgebra.jl/pull/1306 \newunicodechar{└}{\\└} diff --git a/note/.gitignore b/note/.gitignore new file mode 100644 index 0000000..8904018 --- /dev/null +++ b/note/.gitignore @@ -0,0 +1,3 @@ +/.quarto/ +**/*.quarto_ipynb +/.vscode/ \ No newline at end of file diff --git a/note/Note_on_mixed_model_calculations.qmd b/note/Note_on_mixed_model_calculations.qmd new file mode 100644 index 0000000..c81f9c5 --- /dev/null +++ b/note/Note_on_mixed_model_calculations.qmd @@ -0,0 +1,278 @@ +--- +title: "Note on mixed model calculations" +author: + - name: Douglas Bates + email: dmbates@gmail.com + orcid: 0000-0001-8316-9503 + affiliation: + - name: University of Wisconsin - Madison + city: Madison + state: WI + url: https://www.wisc.edu + department: Statistics + - name: Phillip Alday + email: me@phillipalday.com + orcid: 0000-0002-9984-5745 + affiliation: + - name: Beacon Biosignals + url: https://beacon.bio +date: last-modified +date-format: iso +toc: true +bibliography: bibliography.bib +number-sections: true +engine: julia +julia: + exeflags: + - -tauto + - --project=@. +format: + html: + toc: true + toc-location: right + embed-resources: true +--- + +## Introduction {#sec-intro} + +This note is to reproduce the simulated two-factor anova model in @Zhou03042019 and to fit it using the [MixedModels.jl](https://github.com/JuliaStats/MixedModels.jl) package. +I would have included comparative results from the [VarianceComponentModels.jl](https://github.com/OpenMendel/VarianceComponentModels.jl) package but I was unable to install it, even on the long-term-support (lts) version of Julia. +(Repairs may be as simple as merging the [CompatHelper PR](https://github.com/OpenMendel/VarianceComponentModels.jl/pull/22) but I didn't check.) + +The model incorporates two random effects factors and their interaction. +In the simulation the number of levels of the grouping factors for the random effects is set at 5 for `a` and `b` so that the number of levels for the interaction term, `ab`, is 25. +The number of replicates at each level of the interaction is the parameter `c`. + +We note that only having five levels each of the grouping factors $\mathbf{a}$ and $\mathbf{b}$ may result in an unstable estimation situation, because it is difficult to estimate a variance from only five distinct pieces of information. +In such cases, as we will see below, the estimates of variance components can converge to zero. + +Convergence to zero for a variance component is not a problem when using [MixedModels.jl](https://github.com/JuliaStats/MixedModels.jl) because it was designed with this possibility in mind. + +## Setting up the simulation {#sec-setup} + +We create a `DataFrame` corresponding to the experimental design with a placeholder response vector then create a `LinearMixedModel` without fitting it then use the `simulate!` function in the `MixedModels` package to simulate the response given values of the parameter. + +Load the packages to be used + +```{julia} +#| label: loadpackages +using Chairmarks # to benchmark function executions +using DataFrames +using MixedModels +using PooledArrays # similar to factor in R +using Random # random number generators +using RCall # call R from Julia +``` + +Create a function to generate a `DataFrame` with a numeric response, `y`, and `PooledArrays`, `a`, `b`, and `ab`. + +```{julia} +#| label: definedesign +#| output: false +""" + design(c; a=5, b=5) + +Return a DataFrame with columns y, a, b, and ab in a two-factor design with interaction +""" +function design(c::Integer; a::Integer=5, b::Integer=5) + signed, compress = true, true # used as named arguments to PooledArray() + av = repeat(string.('a':'z')[1:a], inner=b * c) + bv = repeat(string.('A':'Z')[1:b], inner=c, outer=b) + return DataFrame( + y = zeros(a * b * c), + a = PooledArray(av; signed, compress), + b = PooledArray(bv; signed, compress), + ab = PooledArray(av .* bv; signed, compress), + ) +end +``` + +For `c = 5` the data frame is + +```{julia} +#| label: d05 +d05 = design(5) +``` + +## Simulating a single instance {#sec-simulating} + +We define a formula, `form`, that provides for scalar fixed-effects and scalar random effects for `a`, `b`, and `ab`. +Because this binding is declared to be `const` we avoid declaring it more than once. + +```{julia} +#| output: false +#| warn: false +if !@isdefined(form) + const form = @formula(y ~ 1 + (1 | a) + (1 | b) + (1 | ab)) +end +m05 = LinearMixedModel(form, d05) +``` + +The `simulate!` function modifies this model by overwriting the place-holder response vector with the simulated response, in the model object only, not in the original data frame, `d05`. +The parameters are set, as in @Zhou03042019, to $\mathbf{\beta}=[1]$, $\sigma=1$ (written $\sigma_e$ in @Zhou03042019) and +values of $\mathbf{\theta}=\left[\sigma_1/\sigma, \sigma_2/\sigma, \sigma_3/\sigma\right]$ corresponding to the (constant)ratios $\sigma^2_i/\sigma^2_e, i=1,2,3$ in @Zhou03042019. + +The first set of $\mathbf{\theta}$ values is zeros, as in @Zhou03042019. + +```{julia} +rng = Xoshiro(6345789) # initialize a random number generator +print(fit!(simulate!(rng, m05; β=[1.0], σ=1.0, θ=zeros(3)))) +``` + +To be able to reproduce this fit we copy the simulated response vector, `m05.y`, into `d05.y`. + +```{julia} +#| label: copym05ytod05y +#| output: false +copyto!(d05.y, m05.y) +``` + +It is not surprising that estimates of two of the three components of $\mathbf{\theta}$ are zero, as the response was simulated from $\mathbf{\theta}=[0,0,0]$. + +The `optsum` property of the fitted model, `m05`, contains information on the progress of the iterative optimization to obtain the mle's. + +```{julia} +#| label: showoptsum +m05.optsum +``` + +In this case it took 55 function evaluations to declare convergence but that number will depend on the convergence criteria set. +We use rather stringent criteria. +In particular, `ftol_rel`, which appears to be the criterion used in @Zhou03042019, is, by default, set to $10^{-12}$ in `MixedModels.jl`, as compared to $10^{-8}$ in @Zhou03042019. + +The progress of the iterations is recorded in the `fitlog` table in the `optsum` property. + +```{julia} +#| label: fitlog +m05.optsum.fitlog +``` + +If we look at the last 15 values of the objective we can see that the optimizer would have been declared to have converged in 15 fewer function evaluations if `ftol_rel` was set to $10^{-8}$. + +```{julia} +#| label: lastobj +last(m05.optsum.fitlog.objective, 15) +``` + +When comparing the number of iterations between algorithms, bear in mind that the optimizer used here reports the number of evaluations of the objective. +Other algorithms may count iterations that involve more than one evaluation of the objective or may involve gradient evaluations. + +## Benchmarking {#sec-benchmark} + +For this model, and for all the other models in the simulation, the blocked Cholesky factor, $\mathbf{L}$ to be updated for each evaluation of the profiled log-likelihood has the structure + +```{julia} +#| label: Blockdesc +BlockDescription(m05) +``` + +That is, regardless of the value of `c` in the simulation, the updates for evaluating the profiled log-likelihood are for a blocked lower Cholesky factor of size $37\times 37$ of which the upper left $25\times 25$ block is diagonal and trivial to update. + +As a sparse matrix this would have the form + +```{julia} +sparseL(m05; full=true) +``` + +but it is stored as a blocked triangular matrix. +The convergence to $\theta_1=0, \theta_2=0$ makes this matrix appear to be more sparse than it, in fact, is. +Generally there would be non-zero values in the last 12 rows and first 30 columns. + +The update operation is very fast + +```{julia} +#| label: bnchmrkupdate +@b m05 objective(updateL!(_)) +``` + +as is the fitting procedure from start to finish + +```{julia} +#| label: bnchmrkfit +@b fit(MixedModel, $form, $d05) +``` + +and the `m05` object is quite small + +```{julia} +#| label: modelsize +Base.summarysize(m05) +``` + +We now fit the same model to this response with the `lme4` package in [R](https://www.r-project.org). + +First, transfer the data frame and the formula to R + +```{julia} +#| output: false +@rput d05 +@rput form +``` + +then fit the model + +```{julia} +#| label: Rfit +#| warning: false +R"m05 <- lme4::lmer(form, d05, REML=FALSE, control=lme4::lmerControl(calc.derivs=FALSE))" +R"summary(m05)" +``` + +Notice that the convergence is to a slightly different parameter vector but a similar deviance. +The main difference in the converged parameter estimates is in $\sigma_2$, which is $1.426\times10^{-5}$ here and zero in the MixedModels fit. +But random effects with a standard deviation that small are negligible. + +```{julia} +#| label: Rfitdeviance +deviance(m05) +``` + +Fitting this model in `R` using `lme4::lmer` is much slower than using `MixedModels.jl`. + +```{julia} +#| warning: false +@b R"lme4::lmer(form, d05, REML=FALSE, lme4::lmerControl(calc.derivs=FALSE))" +``` + +### Data simulated from non-zero $\theta$ values + +```{julia} +fit!(simulate!(rng, m05; β=[1.], σ=1., θ=ones(3))) +``` + +```{julia} +m05.optsum +``` + +## Larger data sets {#sec-larger} + +Increasing `c`, the number of replicates at each level of `ab` does not substantially increase the size of the model + +```{julia} +#| label: m50 +m50 = fit!( + simulate!(rng, LinearMixedModel(form, design(50)); β=ones(1), σ=1.0, θ=ones(3)) +) +print(m50) +``` + +```{julia} +Base.summarysize(m50) +``` + +because the size of `L` and `A` are the same as before + +```{julia} +BlockDescription(m50) +``` + +and the elapsed time per evaluation of the objective is essentially the same as for the smaller model. + +```{julia} +@b m50 objective(updateL!(_)) +``` + +### References {.unnumbered} + +::: {#refs} +::: diff --git a/note/_quarto.yml b/note/_quarto.yml new file mode 100644 index 0000000..475a409 --- /dev/null +++ b/note/_quarto.yml @@ -0,0 +1,3 @@ +project: + title: "Note_on_mixed_model_calculations" + diff --git a/note/bibliography.bib b/note/bibliography.bib new file mode 100644 index 0000000..85609ff --- /dev/null +++ b/note/bibliography.bib @@ -0,0 +1,16 @@ +@article{Zhou03042019, + author = {Hua Zhou and Liuyi Hu and Jin Zhou and Kenneth Lange}, + title = {MM Algorithms for Variance Components Models}, + journal = {Journal of Computational and Graphical Statistics}, + volume = {28}, + number = {2}, + pages = {350--361}, + year = {2019}, + publisher = {ASA Website}, + doi = {10.1080/10618600.2018.1529601}, + note ={PMID: 31592195}, + URL = {https://doi.org/10.1080/10618600.2018.1529601}, + eprint = {https://doi.org/10.1080/10618600.2018.1529601} +} + +