An implementation for SCopt (through Mosek) in Julia

msk_scopt.jl

# This file contains the implementation of the geometric programming interface

export
  MSK_OPR_ENT,
  MSK_OPR_EXP,
  MSK_OPR_LOG,
  MSK_OPR_POW,
  scbegin,
  scend

# A macro to make calling C API a little cleaner
macro scopt_ccall(func, args...)
  f = Base.Meta.quot(symbol("MSK_$(func)"))
  args = [esc(a) for a in args]
  quote
    ccall(($f,libmosekscopt), $(args...))
  end
end

# msk_enums
const MSK_OPR_ENT = convert(Int32,23)
const MSK_OPR_EXP = convert(Int32,1)
const MSK_OPR_LOG = convert(Int32,2)
const MSK_OPR_POW = convert(Int32,3)

type MSKnlhandt
  numcon::Int32
  numvar::Int32
  numopro::Int32
  opro::Vector{Int32}
  oprjo::Vector{Int32}
  oprfo::Vector{Float64}
  oprgo::Vector{Float64}
  oprho::Vector{Float64}

  numoprc::Int32
  oprc::Vector{Int32}
  opric::Vector{Int32}
  oprjc::Vector{Int32}
  oprfc::Vector{Float64}
  oprgc::Vector{Float64}
  oprhc::Vector{Float64}

  ptrc::Vector{Int32}
  subc::Vector{Int32}

  # Work storage employed when evaluating the functions
  ibuf::Vector{Int32}
  zibuf::Vector{Int32}
  zdbuf::Vector{Float64}

  log::Ptr{Void}

  # Error data
  msg::Ptr{Cchar}
  res::Ptr{Int32}
end

# msk_functions
freetask(task::MSKtask, buffer) = @msk_ccall("freetask", Nothing, (Ptr{Void},Ptr{Void}),task,buffer)

function scbegin(task::MSKtask,
                 opro::Vector{Int32},     # Defines the functions used for the objective.
                 oprjo::Vector{Int32},    # Defines the variable indexes used in non-linear objective function.
                 oprfo::Vector{Float64},  # Defines constants used in the objective.
                 oprgo::Vector{Float64},  # Defines constants used in the objective.
                 oprho::Vector{Float64},  # Defines constants used in the objective.
                 oprc::Vector{Int32},     # Defines the functions used for the constraints.
                 opric::Vector{Int32},    # Defines the variable indexes used in the non-linear constraint functions.
                 oprjc::Vector{Int32},    # Defines the constraint indexes where non-linear functions appear.
                 oprfc::Vector{Float64},  # Defines constants used in the non-linear constraints.
                 oprgc::Vector{Float64},  # Defines constants used in the non-linear constraints.
                 oprhc::Vector{Float64})  # Defines constants used in the non-linear constraints.
  res = MSK_RES_ERR_SPACE
  numopro = length(opro)
  numoprc = length(oprc)
  numcon = getnumcon(task)
  numvar = getnumvar(task)
  try
  nlh = MSKnlhandt(numcon,
                   numvar,
                   numopro,
                   Array(Int32, numopro),
                   Array(Int32, numopro),
                   Array(Float64, numopro),
                   Array(Float64, numopro),
                   Array(Float64, numopro),
                   numoprc,
                   Array(Int32, numoprc),
                   Array(Int32, numoprc),
                   Array(Int32, numoprc),
                   Array(Float64, numoprc),
                   Array(Float64, numoprc),
                   Array(Float64, numoprc))
  if numopro > 0
    nlh.opro[:] = opro
    nlh.oprjo[:] = oprjo
    nlh.oprfo[:] = oprfo
    nlh.oprgo[:] = oprgo
    nlh.oprho[:] = oprho
  end
  if numoprc > 0
    nlh.oprc[:] = oprc
    nlh.opric[:] = opric
    nlh.oprjc[:] = oprjc
    nlh.oprfc[:] = oprfc
    nlh.oprgc[:] = oprgc
    nlh.oprhc[:] = oprhc
  end
  nlh.ibuf = Array(Int32, numvar)
  nlh.zibuf = Array(Int32, numvar)
  nlh.zdbuf = Array(Float64, numvar)
  if numoprc > 0
    nlh.ptrc = Array(Int32, numcon+1)
    nlh.subc = Array(Int32, numoprc)
  end
  try
  if (numvar == 0 || nlh.ibuf != C_NULL ) && numoprc == 0
    if numcon > 0 && numoprc > 0
      # Setup of ptrc and sub. Afte the setup then 
      # 1) ptrc[i+1]-ptrc[i]: Is the number of nonlinear terms in the ith constraint.
      # 2) subc[ptrc[i],...,ptrc[i+1]-1]: List the nonlinear terms in the ith constraint.
      for k=1:numoprc
        nlh.ptrc[opric[k]] += 1
      end
      sum = 0
      for k=1:numcon
        nlh.ptrc[k], sum = sum, sum + nlh.ptrc[k]
      end
      for k=1:numoprc
        nlh.subc[nlh.ptrc[opric[k]]] = k
        nlh.ptrc[opric[k]] += 1
      end
      k = nlh.numcon+1
      while k > 1
        nlh.ptrc[k] = nlh.ptrc[k-1]
        k -= 1
      end
      nlh.ptrc[1] = 0
    end
    res = @msk_ccall("putnlfunc",Int32,(Ptr{Void},Ptr{Void},Ptr{Void},
                     Ptr{Void}),task,nlh,scstruc,sceval)
  end
  if res != MSK_RES_OK
    msg = getlasterror(task_)
    throw (MosekError(res,msg))
  end
  sch
end

function scend(task::MSKtask, sch::MSKnlhandt)
  res = @msk_ccall("putnlfunc",Int32,(Ptr{Void},Ptr{Void},Ptr{Void},
                   Ptr{Void}),task,C_NULL,C_NULL,C_NULL)
  if res != MSK_RES_OK
    msg = getlasterror(task_)
    throw (MosekError(res,msg))
  end
end

# The Madness Begins!

# Purpose: Evaluates an operator and its derivatives.
#   fxj_:    Is the function value
#   grdfxj_: Is the first derivative.
#   hexfxj_: Is the second derivative.
function evalopr(opr::Int32,
                 f::Float64,
                 g::Float64,
                 h::Float64,
                 xj::Float64,
                 fxj_::Ptr{Float64},
                 grdfxj_::Ptr{Float64},
                 hesfxj_::Ptr{Float64})
  if opr == MSK_OPR_ENT
    if xj <= 0.0
      return int32(1)::Int32
    end
    if fxj_ != C_NULL
      unsafe_store!(fxj_, f*xj*log(xj))
    end
    if grdfxj_ != C_NULL
      unsafe_store!(grdfxj_, f*(1.0+log(xj)))
    end
    if hesfxj_ != C_NULL
      unsafe_store!(hesfxj_, f/xj)
    end
  elseif opr == MSK_OPR_EXP
    rtemp = exp(g*xj+h)
    if fxj_ != C_NULL
      unsafe_store!(fxj_, f*rtemp)
    end
    if grdfxj_ != C_NULL
      unsafe_store!(grdfxj_, f*g*rtemp)
    end
    if hesfxj_ != C_NULL
      unsafe_store!(hesfxj_, f*g*g*rtemp)
    end
  elseif opr == MSK_OPR_LOG
    rtemp = exp(g*xj+h)
    if rtemp <= 0.0
      return int32(1)::Int32
    end
    if fxj_ != C_NULL
      unsafe_store!(fxj_, f*log(rtemp))
    end
    if grdfxj_ != C_NULL
      unsafe_store!(grdfxj_, (f*g)/rtemp)
    end
    if hesfxj_ != C_NULL
      unsafe_store!(hesfxj_, (f*g*g)/(rtemp^2))
    end
  elseif opr == MSK_OPR_POW
    if fxj_ != C_NULL
      unsafe_store!(fxj_, f*(xj+h)^g)
    end
    if grdfxj_ != C_NULL
      unsafe_store!(grdfxj_, f*g*(xj+h)^(g-1.0))
    end
    if hesfxj_ != C_NULL
      unsafe_store!(hesfxj_, f*g*(g-1.0)*(xj+h)^(g-2.0))
    end
  end
  MSK_RES_OK::Int32
end

# Purpose: Compute number objective value and the gradient.
function scobjeval(nlh::MSKnlhandt,
                   x::Ptr{Float64},
                   objval_::Ptr{Float64},
                   grdnz_::Ptr{Int32},
                   grdsub_::Ptr{Int32},
                   grdval_::Ptr{Float64})
  # temporary objects
  fx = Array(Float64, 1)
  grdfx = Array(Float64, 1)
  objval = Array(Float64, 1)
  grdnz = Array(Float64, 1)
  r = int32(0)

  if objval_ != C_NULL
    objval = pointer_to_array(objval_, 1)
    objval[1] = 0.0
  end
  if grdnz_ != C_NULL
    grdnz = pointer_to_array(grdnz_, 1)
    grdnz[1] = int32(0)
  end
  for k=1:nlh.numopro
    if r != MSK_RES_OK break end
    j = nlh.oprjo[k]
    r = evalopr(nlh.opro[k],nlh.oprfo[k],nlh.oprgo[k],nlh.oprho[k],
                x[j],pointer(fx),pointer(grdfx),C_NULL)
    #should check r == MSK_RES_OK here
    if objval_ != C_NULL
      objval[1] += fx[1]
    end
    if grdnz_ != C_NULL
      nlh.zdbuf[j] += grdfx[1]
      if nlh.zibuf[j] == int32(0)
        # A new nonzero in the gradient has been located
        unsafe_store!(grdsub_, j, grdnz[1])
        nlh.zibuf[j] = int32(1)
        grdnz[1] += 1
      end
    end
  end
  if grdnz != C_NULL
    # Buffers should be zeroed.
    for k=1:grdnz[1]
      j = unsafe_load(grdsub_, k)
      if grdval_ != C_NULL
        unsafe_store!(grdval_, nlh.zdbuf[j], k)
      end
      nlh.zibuf[j] = 0
      nlh.zdbuf[j] = 0.0
    end
  end
  r
end

# Purpose: Compute number value and the gradient of constraint i.
#          grdsub[0,...,grdnz-1] tells which values are needed in gradient
#          that is required.
function scgrdconeval(nlh::MSKnlhandt,
                      i::Int32,
                      x::Ptr{Float64},
                      fval_::Ptr{Float64},
                      grdnz::Int32,
                      grdsub_::Ptr{Int32},
                      grdval_::Ptr{Float64})
  # temporary objects
  fx = Array(Float64, 1)
  grdfx = Array(Float64, 1)
  fval = Array(Float64, 1)
  r = int32(0)
  if fval_ != C_NULL
    fval = pointer_to_array(fval_, 1)
    fval[1] = 0.0
  end
  if nlh.ptrc != C_NULL
    gnz = 0
    for p=nlh.ptrc[i]:nlh.ptrc[i+1]
      if r == MSK_RES_OK break end
      k = nlh.subc[p]
      j = nlh.oprjc[k]
      r = evalopr(nlh.oprc[k],nlh.oprfc[k],nlh.oprgc[k],nlh.oprhc[k],
                  x[j],pointer(fx),pointer(grdfx),C_NULL)
      # should check r == MSK_RES_OK here
      if fval_ != C_NULL
        fval[1] += fx[1]
      end
      if grdnz > 0
        nlh.zdbuf[j] += grdfx[1]
        if nlh.zibuf[j] != 0
          # A new nonzero in the gradient has been located
          nlh.ibuf[gnz] = j
          nlh.zibuf[j] = 1
          gnz += 1
        end
      end
    end
    if grdval_ != C_NULL
      # Setup gradient
      for k=1:grdnz
        j = unsafe_load(grdsub_, k)
        unsafe_store!(grdval_, nlh.zdbuf[j], k)
      end
    end
    for k=1:gnz
      j = nlh.ibuf[k]
      nlh.zibuf[j] = 0
      nlh.zdbuf[j] = 0.0
    end
  elseif grdval_ != C_NULL
    grdval = pointer_to_array(grdval_, grdnz)
    grdval[1:grdnz] = 0.0
  end
  r
end


# Purpose: Evaluate the nonlinear function and return the requested information to MOSEK. 
#  (e.g. function values and in addition 1st and 2nd order derivative information.)
function sceval_wrapper(nlhandle_::Ptr{Void},
                        xx_::Ptr{Float64},
                        yo::Float64,
                        yc_::Ptr{Float64},
                        objval_::Ptr{Float64},
                        numgrdobjnz_::Ptr{Int32},
                        grdobjsub_::Ptr{Int32},
                        grdobjval_::Ptr{Float64},
                        numi::Int32,
                        subi_::Ptr{Int32},
                        conval_::Ptr{Float64},
                        grdconptrb_::Ptr{Int32},
                        grdconptre_::Ptr{Int32},
                        grdconsub_::Ptr{Int32},
                        grdconval_::Ptr{Float64},
                        grdlag_::Ptr{Float64},
                        maxnumhesnz::Int32,
                        numhesnz_::Ptr{Int32},
                        hessubi_::Ptr{Int32},
                        hessubj_::Ptr{Int32},
                        hesval_::Ptr{Float64})
  fx = Array(Float64, 1)
  grdfx = Array(Float64, 1)
  hesfx = Array(Float64, 1)

  nlh = unsafe_pointer_to_objref(nlhandle_)::MSKnlhandt
  numcon = nlh.numcon
  numvar = nlh.numvar
  xx = pointer_to_array(xx_, numvar)
  yc = pointer_to_array(yc_, numcon)

  subi = pointer_to_array(subi_, numi) .+ 1
  grdconptrb = pointer_to_array(grdconptrb_, numi)
  grdconptre = pointer_to_array(grdconptre_, numi)
  numhesnz = pointer_to_array(numhesnz_, 1)

  r = scobjeval(nlh,xx_,objval_,numgrdobjnz_,grdobjsub_,grdobjval_)
  for k=1:numi
    if r != MSK_RES_OK break end
    i = subi[k]
    r = scgrdconeval(nlh,i,xx_,
                     conval_==C_NULL    ? C_NULL : conval_ + (k-1)*sizeof(Float64),
                     grdconsub_==C_NULL ? 0      : grdconptre[k] - grdconptrb[k],
                     grdconsub_==C_NULL ? C_NULL : grdconsub_ + grdconptrb[k]*sizeof(Int32),
                     grdconval_==C_NULL ? C_NULL : grdconval_ + grdconptrb[k]*sizeof(Float64))
  end

  if grdlag_ != C_NULL && r != MSK_RES_OK
    # Compute and store the gradient of the Lagrangian (Note: it is stored as a dense vector)
    grdlag = pointer_to_array(grdlag_, numvar)
    grdlag[1:numvar] = 0.0

    if !isapprox(yo, 0.0)
      for k=1:nlh.numopro
        if r == MSK_RES_OK break end
        j = nlh.oprjo[k];
        r = evalopr(nlh.opro[k],nlh.oprfo[k],nlh.oprgo[k],nlh.oprho[k],
                    xx[j],C_NULL,pointer(grdfx),C_NULL)
        # should check r == MSK_RES_OK here
        grdlag[j] += yo * grdfx[1]
      end
    end
    if nlh.ptrc != C_NULL
      for l=1:numi
        if r != MSK_RES_OK break end
        i = subi[l]
        for p=nlh.ptrc[i]:nlh.ptrc[i+1]
          if r != MSK_RES_OK break end
          k = nlh.subc[p]
          j = nlh.oprjc[k]
          r = evalopr(nlh.oprc[k],nlh.oprfc[k],nlh.oprgc[k],nlh.oprhc[k],
                      xx[j],C_NULL,pointer(grdfx),C_NULL)
          grdlag[j] -= yc[i] * grdfx[1]
        end
      end
    end
  end
  if maxnumhesnz > 0 && r == MSK_RES_OK
    # Compute and store the Hessian of the Lagrange function.
    numhesnz[1] = int32(0)
    if !isapprox(yo_, 0.0)
      for k=1:nlh.numopro
        if r != MSK_RES_OK break end
        j = nlh.oprjo[k]
        r = evalopr(nlh.opro[k],nlh.oprfo[k],nlh.oprgo[k],nlh.oprho[k],
                    xx[j],C_NULL,C_NULL,pointer(hesfx))
        if nlh.zibuf[j] != 0
          numhesnz[1] += 1
          nlh.zibuf[j] = numhesnz[1]
          unsafe_store!(hessubi_, j, nlh.zibuf[j]-1)
          unsafe_store!(hesval_, 0.0, nlh.zibuf[j]-1)
        end
        tmp = unsafe_load(hesval_, zibuf[j]-1)
        unsafe_store!(hesval_, tmp + yo*hesfx[1], zibuf[j]-1)
      end
    end

    if nlh.ptrc != C_NULL
      for l=1:numi
        if r != MSK_RES_OK break end
        i = subi[l]
        for p=nlh.ptrc[i]:nlh.ptrc[i+1]
          if r != MSK_RES_OK break end
          k = nlh.subc[p]
          j = nlh.oprjc[k]
          r = evalopr(nlh.oprc[k],nlh.oprfc[k],nlh.oprgc[k],nlh.oprhc[k],
                      xx[j],C_NULL,C_NULL,pointer(hesfx))
          if nlh.zibuf[j] != 0
            numhesnz[1] += 1
            zibuf[j] = numhesnz[1]
            unsafe_store!(hessubi_, j, nlh.zibuf[j]-1)
            unsafe_store!(hesval_, 0.0, nlh.zibuf[j]-1)
          end
          tmp = unsafe_load(hesval_, zibuf[j]-1)
          unsafe_store!(hesval_, tmp - yc[i]*hesfx[1], zibuf[j]-1)
        end
      end
    end

    if numhesnz[1] > maxnumhesnz
      throw MosekError(r,"Hessian evaluation error")
    end

    for k=1:numhesnz[1]
      j = unsafe_load(hessubi_, k)
      unsafe_store!(hessubj_, j, k)
      nlh.zibuf[j] = 0
    end
  end
  r
end


# Purpose: Provide information to MOSEK about the problem structure and sparsity.
function scstruc_wrapper(nlhandle_::Ptr{Void},
                         numgrdobjnz_::Ptr{Int32},
                         grdobjsub_::Ptr{Int32},
                         i::Int32,
                         convali_::Ptr{Int32},
                         grdconinz_::Ptr{Int32},
                         grdconisub_::Ptr{Int32},
                         yo::Int32,
                         numycnz::Int32,
                         ycsub_::Ptr{Int32},
                         maxnumhesnz::Int32,
                         numhesnz_::Ptr{Int32},
                         hessubi_::Ptr{Int32},
                         hessubj_::Ptr{Int32})
  itemp = Array(Int32, 1)
  nlh = unsafe_pointer_to_objref(nlhandle_)::MSKnlhandt
  if numgrdobjnz_ != C_NULL
    scgrdobjstruc(nlh, numgrdobjnz_, grdobjsub_)
  end
  if convali_ != C_NULL || grdconinz_ != C_NULL
    scgrdconistruc(nlh,i,pointer(itemp),grdconisub_)
    if convali_ != C_NULL
      unsafe_store!(convali_, itemp[1] > 0)
    end
    if grdconinz_ != C_NULL
      unsafe_store!(convali_, itemp[1])
    end
  end
  if numhesnz_ != C_NULL
    numhesnz = pointer_to_array(numhesnz_, 1)
    schesstruc(nlh,yo,numycnz,ycsub,numhesnz,hessubi_)
    # if numhesnz[1] > maxnumhesnz && hessubi_ != C_NULL
    #   # Hessian Size Error
    # end
    if hessubi_ != C_NULL
      # In this case the Hessian is diagonal matrix.
      hessubi = pointer_to_array(hessubi_, numhesnz[1])
      hessubj = pointer_to_array(hessubj_, numhesnz[1])
      hessubj[:] = hessubi
    end
  end
  MSK_RES_OK::Int32
end

# Purpose: Compute number of nonzeros and sparsity pattern
#          of the gradient of the objective function.
function scgrdobjstruc(nlh::MSKnlhandt,
                       nz_::Ptr{Int32},
                       sub_::Ptr{Int32})
  if nz_ != C_NULL
    nz = pointer_to_array(nz_, 1)
    nz[1] = 0
    for k=1:nlh.numopro
      j = nlh.oprjo[k]
      if nlh.zibuf[j] == 0
        # A new nonzero in the gradient of the objective has been located
        if sub_ != C_NULL
          unsafe_store!(sub_, j, nz[1])
        end
        nz[1] += 1
        nlh.zibuf[j] = 1
      end
    end
    # Zero zibuf again
    for k=1:nlh.numopro
      j = nlh.oprjo[k]
      zibuf[j] = 0
    end
  end
end

function scgrdconistruc(nlh::MSKnlhandt,
                        i::Int32,
                        nz_::Ptr{Int32},
                        sub_::Ptr{Int32})
  nz = pointer_to_array(nz_, 1)
  nz[1] = 0
  if pointer(nlh.ptrc) != C_NULL # doublecheck!
    for k=nlh.ptrc[i]:nlh.ptrc[i+1]
      j = nlh.oprjc[nlh.subc[k]]
      if nlh.zibuf[j] == 0
        # A new nonzero in the gradient of the ith constraint has been located
        if sub_ != C_NULL
          unsafe_store!(sub_, j, nz[1])
        end
        nz[1] += 1
        nlh.zibuf[j] = 1
      end
    end
    # Zero zibuf again
    for k=nlh.ptrc[i]:nlh.ptrc[i+1]
      j = nlh.oprjc[nlh.subc[k]]
      zibuf[j] = 0
    end
  end
end


 # Computes the number nonzeros the lower triangular part of
 # the Hessian of the Lagrange function and sparsity pattern. 
 #  nz: Number of nonzeros in the Hessian of the Lagrange function.
 #  sub: List of nonzero diagonal elements in the Hessian.
 #       The separable structure is exploited.
function schesstruc(nlh::MSKnlhandt,
                    yo::Int32,
                    numycnz::Int32,
                    ycsub_::Ptr{Int32},
                    nz_::Ptr{Int32},
                    sub_::Ptr{Int32})
  nz = pointer_to_array(nz_, 1)
  ycsub = pointer_to_array(ycsub_, numycnz)
  nz[1] = 0
  if yo > 0
    # Information about the objective function should be computed
    for k=1:nlh.numopro
      j = nlh.oprjo[k]
      if nlh.zibuf[j] == 0
        # A new nonzero in the gradient has been located
        if sub_ != C_NULL
          unsafe_store!(sub_, j, nz[1])
        end
        nz[1] += 1
        nlh.zibuf[j] = 1
      end
    end
  end
  if nlh.ptrc != C_NULL # doublecheck!
    # Evaluate the sparsity of the Hessian. Only constraints specified
    # by ycsub should be included
    for p=1:numycnz
      i = ycsub[p] # constraint index
      for k=nlh.ptrc[i]:nlh.ptrc[i+1]
        j = nlh.oprjc[nlh.subc[k]]
        if nlh.zibuf[j] == 0
          # A new nonzero diagonal element in the Hessian has been located
          if sub_ != C_NULL
            unsafe_store!(sub_, j, nz[1])
          end
          nz[1] += 1
          nlh.zibuf[j] = 1
        end
      end
    end
  end
  # zero work vectors
  if yo != 0
    for k=1:nlh.numopro
      j = nlh.oprjo[k]
      nlh.zibuf[j] = 0
    end
  end
  if nlh.ptrc != C_NULL # doublecheck!
    for p=1:numycnz
      i = ycsub[p]
      for k=nlh.ptrc[i]:nlh.ptrc[i+1]
        j = nlh.oprjc[nlh.subc[k]]
        nlh.zibuf[j] = 0
      end
    end
  end
  MSK_RES_OK::Int32
end