Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
An implementation for SCopt (through Mosek) in Julia
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.