R/translator.R

In quadmod: Quadratic programming modeling language

make.ids <- structure(function
### Construct matrices and vectors with their unique ids
(...
### Named lists of integer vectors of length <= 2. Each item
### represents an optimization variable and its dimension.
 ){
  vars <- list(...)
  ids <- list()
  i <- 1
  for(N in names(vars)){
    d <- vars[[N]]
    n <- prod(d)
    end <- i+n-1
    a <- structure(array(i:end,d),class="id",variable=N)
    i <- end+1
    ids[[N]] <- a
  }
  ids
### Named list of optimization variable ids, which will be used to
### construct the standard form matrices.
},ex=function(){
  vars <- make.ids(xi=5,beta0=1,beta=2)
  with(vars, xi >= 0)
})

">=.id" <- function(a,constant){
### makes a list of inequality constraints, one for each variable
  a[] >= constant
}

"*.id" <- function(a,constant){
### multiply coefficients, returning a linear combination.
  sum(a[]) * constant
}

"*.lincomb" <- function(vars,x){
### multiply coefficients, returning a linear combination
  if(length(x)==1)x <- rep(x,length(vars))
  stopifnot(length(vars)==length(x))
  for(i in seq_along(vars)){
    vars[[i]]$coef <- vars[[i]]$coef * x[i]
  }
  structure(vars,class="lincomb")
}

"*.lincomblist" <- function(vars,x){
### multiply coefficients, returning a linear combination list
  stopifnot(is.numeric(x))
  if(length(x)==1)x <- rep(x,length(vars))
  stopifnot(length(x)==length(vars))
  for(i in seq_along(vars)){
    vars[[i]] <- vars[[i]] * x[i]
  }
  vars
}

sum.id <- function(a,...){
### sum all variables, returning a linear combination
  sum(a[])
}

">=.lincomb" <- function(vars,constant){
### make an inequality constraint from a linear combination
  arglist <- c(list("geq",constant),vars)
  do.call(constrain,arglist)
}

"==.lincomb" <- function(vars,constant){
### make an equality constraint from a linear combination
  arglist <- c(list("eq",constant),vars)
  do.call(constrain,arglist)
}

">=.lincomblist" <- function(lcs,constant){
### make a list of inequality constraints from a list of linear combinations
  L <- lapply(lcs,function(lc)lc>=constant)
  if(length(L)==1)L <- L[[1]]
  L
}

sum.lincomblist <- function(lcs,...){
### collapse a list of linear combinations into a single linear combination
  structure(lapply(lcs,"[[",1),class="lincomb")
}

"+.lincomb" <- function(v1,v2){
### combine optimization variables in a linear combination
  if(class(v2)%in%c("id","lincomblist"))v2 <- sum(v2)
  structure(c(v1,v2),class="lincomb")
}

### Model: define a less than operator for a linear combination of
### optimization variables, which returns a constraint. Define an
### addition operator for a linear combination, which returns another
### linear combination, which is just a list of variables and
### coefficients. 

"[.id" <- function(var,i,j){
### Make some variables from an id matrix.
  N <- attr(var,"variable")
  linear.combinations <- list()
  if(missing(i))i <- 1:nrow(var)
  if(length(dim(var))==1){
    for(ii in i){
      lc <- structure(list(variable(N,ii)),class="lincomb")
      linear.combinations <- c(linear.combinations,list(lc))
    }
  }else{
    if(missing(j))j <- 1:ncol(var)
    for(ii in i)for(jj in j){
      lc <- structure(list(variable(N,ii,jj)),class="lincomb")
      linear.combinations <- c(linear.combinations,list(lc))
    }
  }
  structure(linear.combinations,class="lincomblist")
}

constrain <- structure(function
### Make a constraint list. Write each constraint as c_1 v_1 + ... +
### c_n v_n (in)equality constant, where c_i are real coefficients and
### v_i are optimization variables.
(equality,
### "eq" or "geq" designating whether this is an equality or an
### inequality constraint.
 constant,
### constant on the right side of the inequality.
 ...
### Variables to include in this constraint, usually using something
### like variable("xi",i,j,-1).
 ){
  variables <- list(...)
  vars <- lapply(variables,structure,class=c("variable","list"))
  structure(list(variables=vars,equality=equality,constant=constant),
            class=c("constraint","list"))
},ex=function(){
  constraints <- list()
  for(i in 1:3){
    constraints <-
      c(constraints,list(constrain("geq",0,variable("xi",i))))
  }

  lapply(1:3,function(i)constrain("geq",0,variable("xi",i)))
})

variable <- structure(function
### Make a list that represents an optimization variable and its
### coefficient in the context of a constraint.
(name,
### Name of the optimization variable, created with make.ids.
 i=1,
### If the optimization variable is a vector, the element. If it is a
### matrix, the row.
 j=NULL,
### Column of the optimization variable.
 coef=1
### Coefficient of the optimization variable.
 ){
  structure(list(variable=name,i=i,j=j,coef=coef),
            class=c("variable","list"))
### List with elements "variable" "i" "j" "coef", of class "variable".
},ex=function(){
  variable("xi",2,coef=-1)
  variable("alpha",3,3)
})

print.variable <- function(x,...){
### print method that shows variable indices and coefficient
  index <- if(is.null(x$j)){
    sprintf("%d",x$i)
  }else{
    sprintf("%d,%d",x$i,x$j)
  }
  cat(sprintf("%10.3f*%s[%s]\n",x$coef,x$variable,index))
}

print.constraint <- function(x,...){
### print method that shows each variable and constant
  for(v in x$variables){
    print(v)
  }
  key <- c(geq=">=",eq="==")
  cat(key[x$equality],x$constant,"\n")
}

is.constraint <- function
### Check if an object is a constraint.
(x
### object to check.
 ){
  "constraint" %in% class(x)
### TRUE if x is of class "constraint".
}

### Reality checks for constraints.
validate.constraint <- function(x,vars){
  problem <- function(msg){
    str(x)
    cat("invalid constraint: ",msg,".\n",sep="")
    stop("constraints did not validate.")
  }
  if(!is.constraint(x))problem('not classed "constraint"')
  if( !(x$equality %in% c("geq","eq")) ){
    problem('equality should be "eq" or "geq"')
  }
  if(!is.numeric(x$constant))problem('constant must be numeric')
  if(length(x$constant) != 1)problem('constant must be length 1')
  if(!is.list(x$variables))problem('variables must be a list')
  v.prob <- function(msg){
    str(vars)
    problem(msg)
  }
  for(v in x$variables){
    if(!is.character(v$variable))v.prob('variable name not character')
    if(length(v$variable) != 1)v.prob('variable name must be length 1')
    if(! v$variable %in% names(vars))v.prob('variable name not in var list')
    ## TODO: verify indices, etc.
  }
}

standard.form.constraints <- structure(function
### Convert constraints to standard form matrix and vector.
(constraints,
### List of constraints created using constrain().
 ids
### List of optimization variable ids created using make.ids().
 ){
  ## check if all constraints are valid wrt ids.
  for(constraint in constraints){
    validate.constraint(constraint,ids)
  }
  ## sort by type of constraint:
  eq <- sapply(constraints,"[[","equality")
  meq <- sum(eq=="eq") ## number of equality constraints
  ordered.constraints <- constraints[order(eq)]
  b0 <- sapply(ordered.constraints,"[[","constant")
  n.constraints <- length(ordered.constraints)
  n.variables <- length(unlist(ids))
  A <- matrix(0,n.variables,n.constraints)
  for(cn in seq_along(ordered.constraints)){
    constraint <- ordered.constraints[[cn]]
    for(v in constraint$variables){
      a <- unclass(ids[[v$variable]])
      if(is.null(a))stop(sprintf("variable %s not found",v$variable))
      A.row <- if(is.null(v$j)){
        a[v$i]
      }else{
        a[v$i,v$j]
      }
      A[A.row,cn] <- v$coef
    }
  }
  list(A=A,b0=b0,meq=meq)
### List of variable coefficients and constraint constants.
},ex=function(){
  ## example: optimization variable \alpha\in\RR^3 subject to the
  ## constraint that it must be on the standard simplex:
  opt.vars <- make.ids(alpha=3)
  constraints <- with(opt.vars,c(alpha >= 0,list(sum(alpha) == 1)))
  standard.form.constraints(constraints,opt.vars)

  ## linear svm example
  set.seed(1)
  p <- 2
  y <- rep(c(-1,1),each=20)
  x <- replicate(p,rnorm(length(y),y))
  plot(x,col=y+2,asp=1)
  n <- nrow(x)
  
  vars <- make.ids(slack=n,intercept=1,normal=p)
  constraints <- vars$slack >= 0
  for(i in 1:n){
    ivars <- with(vars,intercept*y[i] + sum(normal)*(x[i,]*y[i]) + slack[i])
    constraints <- c(constraints,list(ivars >= 1))
  }
  solver.args <- standard.form.constraints(constraints,vars)
  n.vars <- length(unlist(vars))
  Dvec <- rep(1e-6,n.vars)
  Dvec[vars$normal] <- 1
  D <- diag(Dvec)
  d <- rep(0,n.vars)
  d[vars$slack] <- -1 ## C == 1

  sol <- quadprog::solve.QP(D,d,solver.args$A,solver.args$b0)
  slack <- sol$solution[vars$slack]
  normal <- sol$solution[vars$normal]
  intercept <- sol$solution[vars$intercept]

  title(paste("A linear Support Vector Machine (SVM):",
              "margin SVs circled, slack drawn in red for other SVs"))
  abline(-intercept/normal[2],-normal[1]/normal[2])
  abline((1-intercept)/normal[2],-normal[1]/normal[2],lty="dotted")
  abline((-1-intercept)/normal[2],-normal[1]/normal[2],lty="dotted")
  f <- function(x)intercept+sum(normal*x)
  yfx <- apply(x,1,f)*y
  on.margin <- abs(yfx-1)<1e-6
  points(x[on.margin,],cex=2,col=y[on.margin]+2)
  i <- yfx<1
  ## these 2 complicated formulas calculate the point on margin where
  ## the data point starts picking up slack
  x1 <- ((y[i]-intercept)*normal[1]-
         normal[1]*normal[2]*x[i,2]+
         normal[2]^2*x[i,1])/(normal[2]^2+normal[1]^2)
  x2 <- (y[i]-intercept)/normal[2]-normal[1]/normal[2]*x1
  segments(x[i,1],x[i,2],x1,x2,col="red")
  l2norm <- function(x)sqrt(sum(x^2))
  rbind(apply(cbind(x1,x2)-x[i,],1,l2norm)*l2norm(normal),slack[i])
})