build A matrix for quantile_ensemble_flex by combining small pieces

[elray1]

in examples I ran with a large number of training set forecasts, this resulted in a substantial time savings

[ryantibs]

@elray1 Apologies, somehow I never saw this! I must have changed GitHub email notifications at the start of my paternity leave, and I need to revisit those settings.

I must saw I'm pretty surprised that anything involved rbind() is faster than pre-allocating memory and filling it up, as I do. I purposely avoided rbind() because it each time it gets applied (in a loop) it copies over the two matrices in question into a new matrix of larger dimension, and then populates it, which is super wasteful.

But clearly I don't understand the way R works, since somehow pre-allocation is itself not the best option. This article seems to suggest that if you just trick R in recasting data frames as lists, pre-allocation can be made blazingly fast. However, I'm not sure how that would jive with the sparse matrices used in this code.

So I think what I'll try (at some point this next week) is this: I'll try to just issue a single call to sparseMatrix() and pass it (huge) vectors i (row index), j (column vector), x (value). To express the (sparse matrix in triplet notation). I think this should just be a short hop away from the code that's already there, and hopefully much faster that the pre-allocation strategy or your rbind() modification. If you get around to it before me :), I would be happy to look at it, but no worries about that---I can also try it myself once I find a few spare cycles.

[ryantibs]

Dropping this link here, for my own future reference, which seems to support the idea of directly calling sparseMatrix() just once (see top voted answer).

[brookslogan]

Hypothesis regarding pre-allocation vs. rbinding: things work differently when dealing with sparse matrices vs. dense matrices. Details:

• It looks like Matrix favors the dgCMatrix format. Adding an entry to an n x p dgCMatrix requires inserting into the middle of a (np)-length vector/array plus adding to some amount of a p-length vector/array; cost depends on the vector/array implementation. Using rbind probably simplifies the logic and reduces the amount of copying involved, but is still awkward: if entries could be in any of the columns, then adding new rows involves adding to many different locations of the (np)-length array (while updating the p-length array is simply summation). This may be faster with a dgRMatrix, as updating the (np)-length array would be concatenation (and the n-length array would be concatenation with a little addition). A single call to sparseMatrix still seems like it could be the best bet for the largest gains in efficiency, but specifying a dgRMatrix [EDIT: combined with rbinding seems the best bet] might be an easier thing to test out.

[ryantibs]

Thank you @brookslogan! Your answer as to why looping over an inserting things into a sparse matrix could be slow makes a lot of sense.

If you have cycles to try out the sparseMatrix() solution this week, then I'd say go for it. I'm completely swamped unfortunately. And if you want help testing it, then we can probably get help as discussed with Katie on slack.

[ryantibs]

Hi @brookslogan, Can you let me know if you're working on this please? I just want to be aware, so I can know whether to ask for help from somebody else. Thanks.

[brookslogan]

Sorry, I need to sort out how github notifications are arranged in my mailboxes. I can take a look at this tomorrow; do we already have a reprex showing where the current code is slow? (Also, I don't know if this investigation really blocks merging this PR.)

[brookslogan]

Looks like most of the formation overhead can be removed in a default call with n=400, p=200, r=7, and 3 tau groups. In the table below, the rbind_Rsparse4 approach doesn't add too much onto the (conversion to Tsparse and) gurobi call itself, Rsparse_gurobi_only. There are a lot of approaches leading to Rsparse4 due to Matrix operations eagerly converting to Csparse and/or Tsparse when undesired in this situation.

+ + + Unit: milliseconds
                expr       min        lq      mean    median        uq
            prealloc 1775.5012 1780.7368 1807.6820 1798.1763 1808.3637
               rbind  544.6211  552.2913  557.1276  554.3662  563.3712
       rbind_Rsparse  353.0511  364.5689  384.2582  367.2996  375.7139
      rbind_Rsparse2  360.3036  370.7677  389.2553  375.9802  382.4810
    cbind_transpose1  720.2931  721.1751  729.0199  725.4385  735.1734
    cbind_transpose2  711.3897  727.8147  735.1530  733.4693  735.6753
      rbind_Tsparse1  361.6058  364.9112  391.3413  372.7105  381.1248
      rbind_Rsparse3  182.9323  187.7461  195.7835  197.9192  199.4950
      rbind_Rsparse4  175.3971  182.2185  188.1685  187.5419  193.0681
 Rsparse_gurobi_only  149.3352  156.5845  162.1915  164.0331  168.2186
       max neval
 1928.2429    10
  574.2195    10
  523.6527    10
  521.7932    10
  746.2331    10
  778.8293    10
  515.3994    10
  207.6690    10
  210.6774    10
  173.0083    10

rbind_Rsparse4 uses an Rsparse format, but builds everything up into lists and only forms the Matrix object at the very end. One approach I haven't tried yet is to do the same for the Tsparse format (Tsparse1 is still using intermediate Matrix objects); this might be easier to understand and build on, but unlikely to be faster.

There just a little bit more work to be done pasting this into the package and running some characterization/golden tests (already prepared), but a PR with the additional time-save should be ready fairly soon.

Benchmarking code:

library("Matrix") # appears to change the behavior of `t` but not `rbind` as of R 4.0.3; also note that R version >=3.2.0 is required for base rbind to work; rBind might be more reliable across versions but is now deprecated and warns
## May need to add @importFrom Matrix t

devtools::load_all()

quantile_ensemble_flex_prealloc = function(qarr, y, tau, weights, tau_groups,
                                           intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                           noncross=TRUE, q0=NULL,
                                           lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                           params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

  # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$A = Matrix(0, nrow=2*N, ncol=P+N, sparse=TRUE)
  model$sense = model$rhs = rep(NA, 2*N)

  # First part of residual constraint
  for (k in 1:r) {
    model$A[(k-1)*n + 1:n, (k-1)*p + 1:p] = tau[k] * qarr[,,k]
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
  }
  model$A[1:N, P + 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    model$A[(r+k-1)*n + 1:n, (k-1)*p + 1:p] = (tau[k]-1) * qarr[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
  }
  model$A[N + 1:N, P + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    B = Matrix(0, nrow=p*(r-num_labels), ncol=P+N, sparse=TRUE)
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          B[count + (k-1)*p + 1:p, (ind[k]-1)*p + 1:p] = -Ipp
          B[count + (k-1)*p + 1:p, (ind[k+1]-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    B = Matrix(0, nrow=r, ncol=P+N, sparse=TRUE)
    for (k in 1:r) B[k, (k-1)*p + 1:p] = vec
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    B = Matrix(0, nrow=n0*(r-1), ncol=N+P, sparse=TRUE)
    for (k in 1:(r-1)) {
      B[(k-1)*n0 + 1:n0, (k-1)*p + 1:p] = -q0[,,k]
      B[(k-1)*n0 + 1:n0, k*p + 1:p] = q0[,,k+1]
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind = function(qarr, y, tau, weights, tau_groups,
                                        intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                        noncross=TRUE, q0=NULL,
                                        lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                        params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # First part of residual constraint
  for (k in 1:r) {
    B = Matrix(0, nrow=n, ncol=P+N, sparse=TRUE)
    B[1:n, (k-1)*p + 1:p] = tau[k] * qarr[,,k]
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    if (k == 1) {
      model$A = B
    }
    else {
      model$A = rbind(model$A, B)
    }
  }
  model$A[1:N, P + 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    B = Matrix(0, nrow=n, ncol=P+N, sparse=TRUE)
    B[1:n, (k-1)*p + 1:p] = (tau[k]-1) * qarr[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    model$A = rbind(model$A, B)
  }
  model$A[N + 1:N, P + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    B = Matrix(0, nrow=p*(r-num_labels), ncol=P+N, sparse=TRUE)
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          B[count + (k-1)*p + 1:p, (ind[k]-1)*p + 1:p] = -Ipp
          B[count + (k-1)*p + 1:p, (ind[k+1]-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    B = Matrix(0, nrow=r, ncol=P+N, sparse=TRUE)
    for (k in 1:r) B[k, (k-1)*p + 1:p] = vec
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    B = Matrix(0, nrow=n0*(r-1), ncol=N+P, sparse=TRUE)
    for (k in 1:(r-1)) {
      B[(k-1)*n0 + 1:n0, (k-1)*p + 1:p] = -q0[,,k]
      B[(k-1)*n0 + 1:n0, k*p + 1:p] = q0[,,k+1]
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind_RsparseMatrix = function(qarr, y, tau, weights, tau_groups,
                                                intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                noncross=TRUE, q0=NULL,
                                                lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # First part of residual constraint
  for (k in 1:r) {
    B = new("dgRMatrix", Dim = c(n, P+N), p = rep(0L, n+1L))
    B[1:n, (k-1)*p + 1:p] = tau[k] * qarr[,,k]
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    if (k == 1) {
      model$A = B
    }
    else {
      model$A = rbind(model$A, B) # XXX auto-converts back to dgCMatrix...
    }
  }
  model$A[1:N, P + 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    B = new("dgRMatrix", Dim = c(n, P+N), p = rep(0L, n+1L))
    B[1:n, (k-1)*p + 1:p] = (tau[k]-1) * qarr[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    model$A = rbind(model$A, B)
  }
  model$A[N + 1:N, P + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    B = new("dgRMatrix", Dim = c(p*(r-num_labels), P+N), p = rep(0L, p*(r-num_labels)+1L))
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          B[count + (k-1)*p + 1:p, (ind[k]-1)*p + 1:p] = -Ipp
          B[count + (k-1)*p + 1:p, (ind[k+1]-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    B = new("dgRMatrix", Dim = c(r, P+N), p = rep(0L, r+1L))
    for (k in 1:r) B[k, (k-1)*p + 1:p] = vec
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    B = new("dgRMatrix", Dim = c(n0*(r-1L), N+P), p = rep(0L, n0*(r-1L)+1L))
    for (k in 1:(r-1)) {
      B[(k-1)*n0 + 1:n0, (k-1)*p + 1:p] = -q0[,,k]
      B[(k-1)*n0 + 1:n0, k*p + 1:p] = q0[,,k+1]
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind_Rsparse2 = function(qarr, y, tau, weights, tau_groups,
                                                 intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                 noncross=TRUE, q0=NULL,
                                                 lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                 params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  model$A = new("dgRMatrix", Dim = c(0L, P+N), p=0L)

  # First part of residual constraint
  for (k in 1:r) {
    B = new("dgRMatrix", Dim = c(n, P+N), p = rep(0L, n+1L))
    B[1:n, (k-1)*p + 1:p] = tau[k] * qarr[,,k] # XXX appears to convert to dgTMatrix...
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    model$A = rbind(model$A, B) # (noteefficiency: this appears to convert back to dgCMatrix, but building in transpose form is slower)
  }
  model$A[1:N, P + 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    B = new("dgRMatrix", Dim = c(n, P+N), p = rep(0L, n+1L))
    B[1:n, (k-1)*p + 1:p] = (tau[k]-1) * qarr[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    model$A = rbind(model$A, B)
  }
  model$A[N + 1:N, P + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    B = new("dgRMatrix", Dim = c(p*(r-num_labels), P+N), p = rep(0L, p*(r-num_labels)+1L))
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          B[count + (k-1)*p + 1:p, (ind[k]-1)*p + 1:p] = -Ipp
          B[count + (k-1)*p + 1:p, (ind[k+1]-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    B = new("dgRMatrix", Dim = c(r, P+N), p = rep(0L, r+1L))
    for (k in 1:r) B[k, (k-1)*p + 1:p] = vec
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    B = new("dgRMatrix", Dim = c(n0*(r-1L), N+P), p = rep(0L, n0*(r-1L)+1L))
    for (k in 1:(r-1)) {
      B[(k-1)*n0 + 1:n0, (k-1)*p + 1:p] = -q0[,,k]
      B[(k-1)*n0 + 1:n0, k*p + 1:p] = q0[,,k+1]
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  .GlobalEnv[["debug.model"]] <- model
  .GlobalEnv[["debug.params"]] <- params

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

rbind2_dgRMatrix_unchecked_unnamed = function(dgRMatrix1, dgRMatrix2) {
  # note that there are still many checks going on inside `new`; it may be more
  # efficient to wait to convert to dgRMatrix at the last minute
  n1 = length(dgRMatrix1@j)
  Dim1 = dgRMatrix1@Dim
  Dim2 = dgRMatrix2@Dim
  Nrow1 = Dim1[[1L]]
  Nrow2 = Dim2[[1L]]
  # Dimnames1 = dgRMatrix1@Dimnames
  # Dimnames2 = dgRMatrix2@Dimnames
  # Rownames1 = Dimnames1[[1L]]
  # Rownames2 = Dimnames2[[1L]]
  # Rownames = if (is.null(Rownames1)) {
  #              if (is.null(Rownames2)) {
  #                NULL
  #              } else {
  #                c(rep("", Nrow1), Rownames2)
  #              }
  #            } else {
  #              if (is.null(Rownames2)) {
  #                c(Rownames1, rep("", Nrow2))
  #              } else {
  #                c(Rownames1, Rownames2)
  #              }
  #            }
  # Colnames = if (is.null(Dimnames1[[2L]])) Dimnames2[[2L]] else Dimnames1[[2L]]
  new(
    "dgRMatrix",
    j=c(dgRMatrix1@j, dgRMatrix2@j),
    p=c(dgRMatrix1@p, n1 + dgRMatrix2@p[-1L]),
    x=c(dgRMatrix1@x, dgRMatrix2@x),
    Dim=c(Nrow1 + Nrow2, Dim1[[2L]])#, # unchecked for ncol compatibility
    # Dimnames=list(Rownames, Colnames)
  )
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind_Rsparse3 = function(qarr, y, tau, weights, tau_groups,
                                                 intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                 noncross=TRUE, q0=NULL,
                                                 lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                 params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  qarr.213 = aperm(qarr, c(2L,1L,3L))
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # First part of residual constraint
  #
  # This row group of A is of the form cbind(B,I), where B is a
  # "(rectangular-block)-diagonal" matrix where each of the r blocks are n x p,
  # and I is an n*r x n*r identity matrix. Build up the `j` slot for an Rsparse
  # representation within a matrix, then convert to a vector:
  #
  # Start with the entries of j for B correct and for I missing:
  BI.j.mat = matrix(c(seq_len(p)-1L,NA_integer_) + rep(seq.int(0L, (r-1L)*p, p), each=n*(p+1L)), p+1L, N)
  # Fill in entries of j for I:
  BI.j.mat[(p+1L),] <- P + seq_len(N)-1L
  # Also use matrix operations to construct the `x` slot:
  B.x.mat = `dim<-`(qarr.213, c(p*n,r)) %*% Diagonal(r, tau)
  BI.x.mat = rbind(`dim<-`(B.x.mat, c(p,n*r)), 1)
  rm(B.x.mat)
  # Prepare Rsparse BI:
  BI = new("dgRMatrix",
           Dim=c(N, N+P),
           j=as.vector(BI.j.mat),
           p=seq.int(0L, N*(p+1L), p+1L),
           x=as.vector(BI.x.mat))
  rm(BI.j.mat, BI.x.mat)
  # Add to model:
  model$A = BI
  model$sense[1:N] = ">="
  model$rhs[1:N] = rep(y,r) * rep(tau, each=n)
  rm(BI)

  # Second part of residual constraint
  BI.j.mat = matrix(c(seq_len(p)-1L,NA_integer_) + rep(seq.int(0L, (r-1L)*p, p), each=n*(p+1L)), p+1L, N)
  BI.j.mat[(p+1L),] <- P + seq_len(N)-1L
  B.x.mat = `dim<-`(qarr.213, c(p*n,r)) %*% Diagonal(r, tau-1)
  BI.x.mat = rbind(`dim<-`(B.x.mat, c(p,n*r)), 1)
  rm(B.x.mat)
  BI = new("dgRMatrix",
           Dim=c(N, N+P),
           j=as.vector(BI.j.mat),
           p=seq.int(0L, N*(p+1L), p+1L),
           x=as.vector(BI.x.mat))
  rm(BI.j.mat, BI.x.mat)
  model$A = rbind2_dgRMatrix_unchecked_unnamed(model$A, BI)
  model$sense[N + 1:N] = ">="
  model$rhs[N + 1:N] = rep(y,r) * rep(tau-1, each=n)
  rm(BI)

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    next_tau_i_in_group = ave(seq_along(tau_groups), tau_groups, FUN=function(group_member_is) c(group_member_is[-1L], NA_integer_))
    lhs_tau_is = seq_along(tau_groups)[!is.na(next_tau_i_in_group)]
    rhs_tau_is = next_tau_i_in_group [!is.na(next_tau_i_in_group)]
    B.nrow = p*(r-num_labels)
    B = new("dgRMatrix",
            Dim = c(B.nrow, P+N),
            j = as.vector(rbind(rep((lhs_tau_is-1L)*p, each=p) + seq_len(p) - 1L,
                                rep((rhs_tau_is-1L)*p, each=p) + seq_len(p) - 1L)),
            p = seq.int(0L, 2L*B.nrow, 2L),
            x = rep(c(-1,1), B.nrow)
            )
    model$A = rbind2_dgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    B.j = seq_len(P) - 1L
    if (intercept) {
      p_nonintercept = p - 1L
      B.j <- B.j[c(FALSE,rep(TRUE,p_nonintercept))] # (recyles flags to omit r intercept indices)
    } else {
      p_nonintercept = p
    }
    B = new("dgRMatrix",
            Dim = c(r, P+N),
            j = B.j,
            p = seq.int(0L, p_nonintercept*r, p_nonintercept),
            x = rep(1, p_nonintercept*r))
    model$A = rbind2_dgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    q0.2.13 = `dim<-`(aperm(q0, c(2L,1L,3L)), c(p, n0*r))
    B = new("dgRMatrix",
            Dim = c(n0*(r-1L), N+P),
            j = rep(seq.int(0L, p*(r-2L), p), each=2L*n0*p) + seq_len(2L*p) - 1L,
            p = seq.int(0L, 2L*n0*p*(r-1L), 2L*p),
            x = as.vector(rbind(-q0.2.13[,-n0*r+seq_len(n0)-1L], q0.2.13[,-seq_len(n0)])))
    model$A = rbind2_dgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

rbind2_listdgRMatrix_unchecked_unnamed = function(listdgRMatrix1, listdgRMatrix2) {
                                        # note that there are still many checks going on inside `new`; it may be more
                                        # efficient to wait to convert to dgRMatrix at the last minute
  n1 = length(listdgRMatrix1[["j"]])
  Dim1 = listdgRMatrix1[["Dim"]]
  Dim2 = listdgRMatrix2[["Dim"]]
  Nrow1 = Dim1[[1L]]
  Nrow2 = Dim2[[1L]]
  list(
    Dim=c(Nrow1 + Nrow2, Dim1[[2L]]), # unchecked for ncol compatibility
    j=c(listdgRMatrix1[["j"]], listdgRMatrix2[["j"]]),
    p=c(listdgRMatrix1[["p"]], n1 + listdgRMatrix2[["p"]][-1L]),
    x=c(listdgRMatrix1[["x"]], listdgRMatrix2[["x"]])
  )
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind_Rsparse4 = function(qarr, y, tau, weights, tau_groups,
                                                 intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                 noncross=TRUE, q0=NULL,
                                                 lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                 params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  qarr.213 = aperm(qarr, c(2L,1L,3L))
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # First part of residual constraint
  #
  # This row group of A is of the form cbind(B,I), where B is a
  # "(rectangular-block)-diagonal" matrix where each of the r blocks are n x p,
  # and I is an n*r x n*r identity matrix. Build up the `j` slot for an Rsparse
  # representation within a matrix, then convert to a vector:
  #
  # Start with the entries of j for B correct and for I missing:
  BI.j.mat = matrix(c(seq_len(p)-1L,NA_integer_) + rep(seq.int(0L, (r-1L)*p, p), each=n*(p+1L)), p+1L, N)
  # Fill in entries of j for I:
  BI.j.mat[(p+1L),] <- P + seq_len(N)-1L
  # Also use matrix operations to construct the `x` slot:
  B.x.mat = `dim<-`(qarr.213, c(p*n,r)) %*% Diagonal(r, tau)
  BI.x.mat = rbind(`dim<-`(B.x.mat, c(p,n*r)), 1)
  rm(B.x.mat)
  # Prepare Rsparse BI:
  BI = list(Dim=c(N, N+P),
            j=as.vector(BI.j.mat),
            p=seq.int(0L, N*(p+1L), p+1L),
            x=as.vector(BI.x.mat))
  rm(BI.j.mat, BI.x.mat)
  # Add to model:
  model$A = BI
  model$sense[1:N] = ">="
  model$rhs[1:N] = rep(y,r) * rep(tau, each=n)
  rm(BI)

  # Second part of residual constraint
  BI.j.mat = matrix(c(seq_len(p)-1L,NA_integer_) + rep(seq.int(0L, (r-1L)*p, p), each=n*(p+1L)), p+1L, N)
  BI.j.mat[(p+1L),] <- P + seq_len(N)-1L
  B.x.mat = `dim<-`(qarr.213, c(p*n,r)) %*% Diagonal(r, tau-1)
  BI.x.mat = rbind(`dim<-`(B.x.mat, c(p,n*r)), 1)
  rm(B.x.mat)
  BI = list(Dim=c(N, N+P),
            j=as.vector(BI.j.mat),
            p=seq.int(0L, N*(p+1L), p+1L),
            x=as.vector(BI.x.mat))
  rm(BI.j.mat, BI.x.mat)
  model$A = rbind2_listdgRMatrix_unchecked_unnamed(model$A, BI)
  model$sense[N + 1:N] = ">="
  model$rhs[N + 1:N] = rep(y,r) * rep(tau-1, each=n)
  rm(BI)

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    next_tau_i_in_group = ave(seq_along(tau_groups), tau_groups, FUN=function(group_member_is) c(group_member_is[-1L], NA_integer_))
    lhs_tau_is = seq_along(tau_groups)[!is.na(next_tau_i_in_group)]
    rhs_tau_is = next_tau_i_in_group [!is.na(next_tau_i_in_group)]
    B.nrow = p*(r-num_labels)
    B = list(Dim = c(B.nrow, P+N),
             j = as.vector(rbind(rep((lhs_tau_is-1L)*p, each=p) + seq_len(p) - 1L,
                                 rep((rhs_tau_is-1L)*p, each=p) + seq_len(p) - 1L)),
             p = seq.int(0L, 2L*B.nrow, 2L),
             x = rep(c(-1,1), B.nrow)
             )
    model$A = rbind2_listdgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    B.j = seq_len(P) - 1L
    if (intercept) {
      p_nonintercept = p - 1L
      B.j <- B.j[c(FALSE,rep(TRUE,p_nonintercept))] # (recyles flags to omit r intercept indices)
    } else {
      p_nonintercept = p
    }
    B = list(Dim = c(r, P+N),
             j = B.j,
             p = seq.int(0L, p_nonintercept*r, p_nonintercept),
             x = rep(1, p_nonintercept*r))
    model$A = rbind2_listdgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    q0.2.13 = `dim<-`(aperm(q0, c(2L,1L,3L)), c(p, n0*r))
    B = list(Dim = c(n0*(r-1L), N+P),
             j = rep(seq.int(0L, p*(r-2L), p), each=2L*n0*p) + seq_len(2L*p) - 1L,
             p = seq.int(0L, 2L*n0*p*(r-1L), 2L*p),
             x = as.vector(rbind(-q0.2.13[,-n0*r+seq_len(n0)-1L], q0.2.13[,-seq_len(n0)])))
    model$A = rbind2_listdgRMatrix_unchecked_unnamed(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Transformations common to both solvers:
  model$A <- new("dgRMatrix", Dim=model$A[["Dim"]], j=model$A[["j"]], p=model$A[["p"]], x=model$A[["x"]])

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_rbind_Tsparse1 = function(qarr, y, tau, weights, tau_groups,
                                                 intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                 noncross=TRUE, q0=NULL,
                                                 lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                 params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  model$A = new("dgTMatrix", Dim = c(0L, P+N))

  # First part of residual constraint
  for (k in 1:r) {
    B = new("dgTMatrix", Dim = c(n, P+N))
    B[1:n, (k-1)*p + 1:p] = tau[k] * qarr[,,k]
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    model$A = rbind(model$A, B) # (noteefficiency: this appears to convert back to dgCMatrix, but building in transpose form is slower)
  }
  model$A[1:N, P + 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    B = new("dgTMatrix", Dim = c(n, P+N))
    B[1:n, (k-1)*p + 1:p] = (tau[k]-1) * qarr[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    model$A = rbind(model$A, B)
  }
  model$A[N + 1:N, P + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    B = new("dgTMatrix", Dim = c(p*(r-num_labels), P+N))
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          B[count + (k-1)*p + 1:p, (ind[k]-1)*p + 1:p] = -Ipp
          B[count + (k-1)*p + 1:p, (ind[k+1]-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    B = new("dgTMatrix", Dim = c(r, P+N))
    for (k in 1:r) B[k, (k-1)*p + 1:p] = vec
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    B = new("dgTMatrix", Dim = c(n0*(r-1L), N+P))
    for (k in 1:(r-1)) {
      B[(k-1)*n0 + 1:n0, (k-1)*p + 1:p] = -q0[,,k]
      B[(k-1)*n0 + 1:n0, k*p + 1:p] = q0[,,k+1]
    }
    model$A = rbind(model$A, B)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_cbind_transpose1 = function(qarr, y, tau, weights, tau_groups,
                                                  intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                  noncross=TRUE, q0=NULL,
                                                  lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                  params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # It is natural to build up the A matrix by `rbind`ing on additional (groups
  # of) rows onto it. However, the Matrix package heavily favors (often converts
  # to) Csparse formats, so instead build up AT, the transpose of A, by
  # `cbind`ing on additional contents, then transpose AT to get A at the end.
  #
  # (noteefficiency: should also try building larger sections the matrix at once, specifying directly in Csparse format, building up in Tsparse format (using Matrix Tsparse objects along the way or only at the end), and/or implementing tau groups with a reduced number of variables and constraints)
  AT = new("dgCMatrix", Dim = c(P+N, 0L), i=integer(0L), x=numeric(0L), p=0L)

  # First part of residual constraint
  for (k in 1:r) {
    BT = new("dgCMatrix", Dim=c(P+N, n), p=rep(0L, n+1L))
    BT[(k-1)*p + 1:p, 1:n] = tau[[k]] * t(qarr[,,k]) # (noteefficiency: should also try aperm'ing qarr at the beginning, as an alternative to multiple `t` calls (examining both computation time and space), and maybe also scaling an aperm'd qarr as a whole rather than in sections (but any potential benefit will probably come from enabling removal of loop overhead (likely minimal) or potentially facilitating more efficient arrangements of A))
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    AT = cbind(AT, BT)
  }
  AT[P + 1:N, 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    BT = new("dgCMatrix", Dim = c(P+N, n), p = rep(0L, n+1L))
    BT[(k-1)*p + 1:p, 1:n] = (tau[[k]]-1) * t(qarr[,,k])
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    AT = cbind(AT, BT)
  }
  AT[P + 1:N, N + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    BT = new("dgCMatrix", Dim = c(P+N, p*(r-num_labels)), p = rep(0L, p*(r-num_labels)+1L))
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          BT[(ind[k]-1)*p + 1:p, count + (k-1)*p + 1:p] = -Ipp
          BT[(ind[k+1]-1)*p + 1:p, count + (k-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    ## B = new("dgRMatrix", Dim = c(r, P+N), p = rep(0L, r+1L))
    BT = new("dgCMatrix", Dim = c(P+N, r), p = rep(0L, r+1L))
    for (k in 1:r) BT[(k-1)*p + 1:p, k] = vec
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    BT = new("dgCMatrix", Dim = c(N+P, n0*(r-1L)), p = rep(0L, n0*(r-1L)+1L))
    for (k in 1:(r-1)) {
      BT[(k-1)*p + 1:p, (k-1)*n0 + 1:n0] = -t(q0[,,k])
      BT[k*p + 1:p, (k-1)*n0 + 1:n0] = t(q0[,,k+1])
    }
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Construct A:
  model$A = t(AT) # (noteefficiency: should also try converting to Tsparse first to see if it's faster)
  rm(AT) # in case AT is large (if so, might want to manually gc here)

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

## building off of approach from @elray1 in https://github.com/ryantibs/quantgen/pull/6
quantile_ensemble_flex_cbind_transpose2 = function(qarr, y, tau, weights, tau_groups,
                                                   intercept=FALSE, nonneg=TRUE, unit_sum=TRUE,
                                                   noncross=TRUE, q0=NULL,
                                                   lp_solver=c("glpk","gurobi"), time_limit=NULL,
                                                   params=list(), verbose=FALSE) {
  # Set up some basic objects that we will need
  n = dim(qarr)[1]
  p = dim(qarr)[2]
  r = dim(qarr)[3]
  N = n*r; P=p*r
  INN = Diagonal(N)
  qarr213 = aperm(qarr, c(2L,1L,3L))
  model = list()

  # Determine LP solver
  use_gurobi = FALSE
  if (lp_solver == "gurobi") {
    if (requireNamespace("gurobi", quietly=TRUE)) use_gurobi = TRUE
    else warning("gurobi R package not installed, using Rglpk instead.")
  }

 # Gurobi setup
  if (use_gurobi) {
    if (!is.null(time_limit)) params$TimeLimit = time_limit
    if (is.null(params$LogToConsole)) params$LogToConsole = 0
    if (verbose) params$LogToConsole = 1
    equal_sign = "="
  }

  # GLPK setup
  else {
    if (!is.null(time_limit)) params$tm_limit = time_limit * 1000
    if (verbose) params$verbose = TRUE
    equal_sign = "=="
  }

  # Vector of objective coefficients
  model$obj = c(rep(0,P), rep(weights, each=r))

  # Matrix of constraint coefficients
  model$sense = model$rhs = rep(NA, 2*N)

  # It is natural to build up the A matrix by `rbind`ing on additional (groups
  # of) rows onto it. However, the Matrix package heavily favors (often converts
  # to) Csparse formats, so instead build up AT, the transpose of A, by
  # `cbind`ing on additional contents, then transpose AT to get A at the end.
  #
  # (noteefficiency: should also try building larger sections the matrix at once, specifying directly in Csparse format, building up in Tsparse format (using Matrix Tsparse objects along the way or only at the end), and/or implementing tau groups with a reduced number of variables and constraints)
  AT = new("dgCMatrix", Dim = c(P+N, 0L), i=integer(0L), x=numeric(0L), p=0L)

  # First part of residual constraint
  for (k in 1:r) {
    BT = new("dgCMatrix", Dim=c(P+N, n), p=rep(0L, n+1L))
    BT[(k-1)*p + 1:p, 1:n] = tau[[k]] * qarr213[,,k] # (noteefficiency: should maybe try different aperm and perhaps scaling as a whole rather than in sections (but any potential benefit will probably come from enabling removal of loop overhead (likely minimal) or potentially facilitating more efficient arrangements of A))
    model$sense[(k-1)*n + 1:n] = ">="
    model$rhs[(k-1)*n + 1:n] = tau[k] * y
    AT = cbind(AT, BT)
  }
  AT[P + 1:N, 1:N] = INN

  # Second part of residual constraint
  for (k in 1:r) {
    BT = new("dgCMatrix", Dim = c(P+N, n), p = rep(0L, n+1L))
    BT[(k-1)*p + 1:p, 1:n] = (tau[[k]]-1) * qarr213[,,k]
    model$sense[(r+k-1)*n + 1:n] = ">="
    model$rhs[(r+k-1)*n + 1:n] = (tau[k]-1) * y
    AT = cbind(AT, BT)
  }
  AT[P + 1:N, N + 1:N] = INN

  # Group equality constraints, if needed
  labels = unique(tau_groups); num_labels = length(labels)
  if (num_labels < r) {
    BT = new("dgCMatrix", Dim = c(P+N, p*(r-num_labels)), p = rep(0L, p*(r-num_labels)+1L))
    Ipp = Diagonal(p)
    count = 0
    for (label in labels) {
      ind = which(tau_groups == label)
      if (length(ind) > 1) {
        for (k in 1:(length(ind)-1)) {
          BT[(ind[k]-1)*p + 1:p, count + (k-1)*p + 1:p] = -Ipp
          BT[(ind[k+1]-1)*p + 1:p, count + (k-1)*p + 1:p] = Ipp
        }
        count = count + (length(ind)-1)*p
      }
    }
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(equal_sign, p*(r-num_labels)))
    model$rhs = c(model$rhs, rep(0, p*(r-num_labels)))
  }

  # Unit sum constraints on alpha, if we're asked to
  if (unit_sum) {
    vec = rep(1,p); if (intercept) vec[1] = 0
    BT = new("dgCMatrix", Dim = c(P+N, r), p = rep(0L, r+1L))
    for (k in 1:r) BT[(k-1)*p + 1:p, k] = vec
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(equal_sign, r))
    model$rhs = c(model$rhs, rep(1, r))
  }

  # Remove nonnegativity constraint on alpha, if we're asked to
  if (!nonneg) {
    if (use_gurobi) model$lb = c(rep(-Inf,P), rep(0,N))
    else model$lb = list(lower=list(ind=1:P, val=rep(-Inf,P)))
  }

  # Remove nonnegativity constraint on intercepts, if needed
  if (intercept && nonneg) {
    if (use_gurobi) model$lb = c(rep(c(-Inf, rep(0,p-1)), r), rep(0,N))
    else model$lb = list(lower=list(ind=(0:(r-1))*p + 1, val=rep(-Inf,r)))
  }

  # Noncrossing constraints, if we're asked to
  if (noncross) {
    n0 = dim(q0)[1]
    q0.213 = aperm(q0, c(2L,1L,3L))
    BT = new("dgCMatrix", Dim = c(N+P, n0*(r-1L)), p = rep(0L, n0*(r-1L)+1L))
    for (k in 1:(r-1)) {
      BT[(k-1)*p + 1:p, (k-1)*n0 + 1:n0] = -q0.213[,,k]
      BT[k*p + 1:p, (k-1)*n0 + 1:n0] = q0.213[,,k+1]
    }
    AT = cbind(AT, BT)
    model$sense = c(model$sense, rep(">=", n0*(r-1)))
    model$rhs = c(model$rhs, rep(0, n0*(r-1)))
  }

  # Construct A:
  model$A = t(AT) # (noteefficiency: should also try converting to Tsparse first to see if it's faster)
  rm(AT) # in case AT is large (if so, might want to manually gc here)

  # Call Gurobi's LP solver, store results
  if (use_gurobi) {
    a = gurobi::gurobi(model=model, params=params)
    alpha = matrix(a$x[1:P],p,r)
    status = a$status
  }

  # Call GLPK's LP solver, store results
  else {
    ## Rglpk seems to recognize Tsparse and Csparse matrices, but not Rsparse; convert beforehand:
    a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
                       rhs=model$rhs, bounds=model$lb, control=params)
    alpha = matrix(a$solution[1:P],p,r)
    status = a$status
  }

  return(enlist(alpha, status))
}

withr::with_rng_version("1.6.0", withr::with_seed(222017L, {
  ## Set up fake problem:
  n = 400L # number of training instances
  p = 20L # number of ensemble components
  tau = c(0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99)
  tau_groups = c(1L, 1L, 2L, 2L, 2L, 3L, 3L)
  ##
  ## n = 3L # number of training instances
  ## p = 2L # number of ensemble components
  ## tau = c(0.1,0.7,0.8,0.9)
  ## tau_groups = c(1L,2L,2L,2L)
  ##
  ## n = 3L # number of training instances
  ## p = 2L # number of ensemble components
  ## tau = c(0.1,0.9)
  ## tau_groups = c(1L,2L)
  ##
  r = length(tau) # number of quantile levels
  qarr = aperm(apply(array(rexp(n * p * r), c(n, p, r)), 1:2, sort), c(2:3, 1L))
  y = rowSums(rnorm(p * r, , 5) * qarr + rnorm(n * p * r) * qarr) + rnorm(n)
}))

quantile_ensemble_prealloc = quantile_ensemble; mockery::stub(quantile_ensemble_prealloc, "quantile_ensemble_flex", quantile_ensemble_flex_prealloc)
quantile_ensemble_rbind = quantile_ensemble; mockery::stub(quantile_ensemble_rbind, "quantile_ensemble_flex", quantile_ensemble_flex_rbind)
quantile_ensemble_rbind_RsparseMatrix = quantile_ensemble; mockery::stub(quantile_ensemble_rbind_RsparseMatrix, "quantile_ensemble_flex", quantile_ensemble_flex_rbind_RsparseMatrix)
quantile_ensemble_rbind_Rsparse2 = quantile_ensemble; mockery::stub(quantile_ensemble_rbind_Rsparse2, "quantile_ensemble_flex", quantile_ensemble_flex_rbind_Rsparse2)
quantile_ensemble_rbind_Rsparse3 = quantile_ensemble; mockery::stub(quantile_ensemble_rbind_Rsparse3, "quantile_ensemble_flex", quantile_ensemble_flex_rbind_Rsparse3)
quantile_ensemble_rbind_Rsparse4 = quantile_ensemble; mockery::stub(quantile_ensemble_rbind_Rsparse4, "quantile_ensemble_flex", quantile_ensemble_flex_rbind_Rsparse4)
quantile_ensemble_rbind_Tsparse1 = quantile_ensemble; mockery::stub(quantile_ensemble_rbind_Tsparse1, "quantile_ensemble_flex", quantile_ensemble_flex_rbind_Tsparse1)
quantile_ensemble_cbind_transpose1 = quantile_ensemble; mockery::stub(quantile_ensemble_cbind_transpose1, "quantile_ensemble_flex", quantile_ensemble_flex_cbind_transpose1)
quantile_ensemble_cbind_transpose2 = quantile_ensemble; mockery::stub(quantile_ensemble_cbind_transpose2, "quantile_ensemble_flex", quantile_ensemble_flex_cbind_transpose2)

result_prealloc = quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups)
result_rbind = quantile_ensemble_rbind(qarr,y,tau,tau_groups=tau_groups)
result_rbind_RsparseMatrix = quantile_ensemble_rbind_RsparseMatrix(qarr,y,tau,tau_groups=tau_groups)
result_rbind_Rsparse2 = quantile_ensemble_rbind_Rsparse2(qarr,y,tau,tau_groups=tau_groups)
result_rbind_Rsparse3 = quantile_ensemble_rbind_Rsparse3(qarr,y,tau,tau_groups=tau_groups)
result_rbind_Rsparse4 = quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups)
result_rbind_Tsparse1 = quantile_ensemble_rbind_Tsparse1(qarr,y,tau,tau_groups=tau_groups)
result_cbind_transpose1 = quantile_ensemble_cbind_transpose1(qarr,y,tau,tau_groups=tau_groups)
result_cbind_transpose2 = quantile_ensemble_cbind_transpose2(qarr,y,tau,tau_groups=tau_groups)

stopifnot(all.equal(result_prealloc, result_rbind))
stopifnot(all.equal(result_prealloc, result_rbind_RsparseMatrix))
stopifnot(all.equal(result_prealloc, result_rbind_Rsparse2))
stopifnot(all.equal(result_prealloc, result_rbind_Rsparse3))
stopifnot(all.equal(result_prealloc, result_rbind_Rsparse4))
stopifnot(all.equal(result_prealloc, result_rbind_Tsparse1))
stopifnot(all.equal(result_prealloc, result_cbind_transpose1))
stopifnot(all.equal(result_prealloc, result_cbind_transpose2))

stopifnot(all.equal(quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,intercept=TRUE,noncross=FALSE),
                    quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,intercept=TRUE,noncross=FALSE)))

stopifnot(all.equal(quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,noncross=FALSE),
                    quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,noncross=FALSE)))

stopifnot(all.equal(quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,unit_sum=FALSE),
                    quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,unit_sum=FALSE)))

result_gurobi_prealloc = quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind = quantile_ensemble_rbind(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind_RsparseMatrix = quantile_ensemble_rbind_RsparseMatrix(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind_Rsparse2 = quantile_ensemble_rbind_Rsparse2(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind_Rsparse3 = quantile_ensemble_rbind_Rsparse3(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind_Rsparse4 = quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_rbind_Tsparse1 = quantile_ensemble_rbind_Tsparse1(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_cbind_transpose1 = quantile_ensemble_cbind_transpose1(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")
result_gurobi_cbind_transpose2 = quantile_ensemble_cbind_transpose2(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi")

stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind_RsparseMatrix))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind_Rsparse2))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind_Rsparse3))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind_Rsparse4))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_rbind_Tsparse1))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_cbind_transpose1))
stopifnot(all.equal(result_gurobi_prealloc, result_gurobi_cbind_transpose2))
## more extensive testing should be performed if an alternative method is selected

microbenchmark::microbenchmark(
                  ## quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind_RsparseMatrix(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind_Rsparse2(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind_Rsparse3(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_rbind_Tsparse1(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_cbind_transpose1(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  ## quantile_ensemble_cbind_transpose2(qarr,y,tau,tau_groups=tau_groups,lp_solver="glpk"),
                  prealloc=quantile_ensemble_prealloc(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind=quantile_ensemble_rbind(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind_Rsparse=quantile_ensemble_rbind_RsparseMatrix(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind_Rsparse2=quantile_ensemble_rbind_Rsparse2(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  cbind_transpose1=quantile_ensemble_cbind_transpose1(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  cbind_transpose2=quantile_ensemble_cbind_transpose2(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind_Tsparse1=quantile_ensemble_rbind_Tsparse1(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind_Rsparse3=quantile_ensemble_rbind_Rsparse3(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  rbind_Rsparse4=quantile_ensemble_rbind_Rsparse4(qarr,y,tau,tau_groups=tau_groups,lp_solver="gurobi"),
                  Rsparse_gurobi_only=gurobi::gurobi(model=debug.model, params=debug.params),
                  times=10L)

[ryantibs]

Thanks @brookslogan, this looks very thorough! Great job. Can't tell easily based on your naming conventions, but were any of these the "just one big call to sparseMatrix()" approach?

Looking forward to seeing wrapped and will look out for your PR. Did you branch off of this this branch? If so then you can actually probably just commit straight to this PR, no? Unless @elray1 wants this merged right away for the ensemble work (which I assumed wasn't the case, based on my other conversations with him, since he's currently moved over to his own Python-based implementations anyway).

[brookslogan]

No, the "one big call to sparseMatrix" would be the same as or one step beyond the currently unimplemented Tsparse-outside-of-Matrix-objects approach. (The extra step would be pre-allocating or building i,j,x vectors of the appropriate size in one go rather than concatenating many times.)

rbind_Rsparse4 is fairly close, doing Rsparse-outside-of-Matrix-objects. (It's not doing that extra pre-allocation/few-concatenations step either.)

I haven't tested whether I can just straight push commits to this PR but I will test it out.