R/runOptimizationOptiSolve.R

In dp3ll1n/SLCDP: Simulation and inference of cell differentiation process governed by logistic growth dynamics

#' Calculate the solution for the Quadratic Programming (QP) problem using Cplex optimization environment.
#'
#' It combines input matrices and constraints for solving the QP problem.
#' @param MX Predictor matrix. It can be M matrix in OLS or Jacobian matrix in the Gauss-Newton algorithm.
#' @param dX Response vector. It can be deltas vector OLS or endValues in the Gauss-Newton algorithm.
#' @param W1Inv Inverse of the estimated covariance matrix. If NULL it is set to the appropriate identity matrix.
#' @param constraints list containing the constraints as generated by setConstraints.
#' @param bvec=NULL values for the inequality constraints (rates >= bvec). If NULL, bvec is filled with 0s.
#' @return  QP solution.

#' @export
#' @examples
#' runOptimizationOptiSolve(MX,dX,W1Inv=NULL,constraints,bvec=NULL)


runOptimizationOptiSolve=function(MX,dX,W1Inv=NULL,constraints,bvec=NULL){
  if(is.null(W1Inv)){
    W1Inv= Matrix::Diagonal(nrow(MX), x = 1)
  }
  if(is.null(bvec)){
    bvec= rep(0, nrow(constraints$A))
  }

  Qmat = (t(MX)%*% W1Inv %*%MX);
  cvec = -2*(t(as.vector(dX)))%*%W1Inv%*%MX;
  cvec=as.matrix(cvec)
  Qmat=as.matrix(Qmat)
  mycop <- optiSolve::cop(f  = quadfun(Q=Matrix::forceSymmetric(Qmat), a=cvec, d=0),
               lb = lbcon(val=constraints$lb),
               ub = ubcon(val=constraints$ub),
               lc = lincon(A= (-1*constraints$A),
                           dir=rep("<=",nrow(constraints$A) ),
                           val=bvec,
                           name=seq_along(1:nrow(constraints$A))))

  res <- optiSolve::solvecop(mycop, solver="alabama", make.definite=T, quiet=FALSE)
  a=res$x
 # a[a<constraints$lb]=constraints$lb[a<constraints$lb]
 # a=zapsmall(a)
  return(a)
}