From a06f60a47bd5f5106b9b4bbb26cf54faf0acaa68 Mon Sep 17 00:00:00 2001 From: Thibaut Lamadon Date: Thu, 14 Mar 2013 15:41:10 +0000 Subject: [PATCH] finished coding the piping between mpeccable and ipopt. I am still getting an error when calling ipopt on the return value of the jacobian. --- DESCRIPTION | 1 + NAMESPACE | 3 + R/optim.variables.r | 8 +- R/solve.mpec.r | 102 +++++++++++++++++++ R/utils.R | 14 +++ examples/mpeccable.csolve.ex.r | 12 +-- test/solve.mpec.r | 175 --------------------------------- 7 files changed, 132 insertions(+), 183 deletions(-) create mode 100644 R/solve.mpec.r delete mode 100644 test/solve.mpec.r diff --git a/DESCRIPTION b/DESCRIPTION index 47351da..5261900 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,3 +29,4 @@ Collate: 'utilsSpline.R' 'SolveDiff.R' 'optim.variables.r' + 'solve.mpec.r' diff --git a/NAMESPACE b/NAMESPACE index 601a30f..2795dee 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,9 +4,12 @@ export(F_Spline2D) export(F_Spline3D) export(F_SplineInt1D) export(computeVarRanges) +export(getSparseValues) export(mpec.vars.collate) export(mpec.vars.desc) export(mpec.vars.extract) +export(mpec.vars.init) +export(mpeccable.csolve) export(param0) exportClasses(FDiff) exportMethods("%*%") diff --git a/R/optim.variables.r b/R/optim.variables.r index 21ddbc4..135ecec 100644 --- a/R/optim.variables.r +++ b/R/optim.variables.r @@ -37,4 +37,10 @@ mpec.vars.collate <- function(x0,vars,coloring=FALSE) { names(ll) = names(vars) return(ll) -} \ No newline at end of file +} + +#' create the list of params, correctly initialized +#' @export +mpec.vars.init <- function(vars,value=0) { + return(rep(value,Reduce('+', vars))) +} diff --git a/R/solve.mpec.r b/R/solve.mpec.r new file mode 100644 index 0000000..3971dbb --- /dev/null +++ b/R/solve.mpec.r @@ -0,0 +1,102 @@ +#' solves a set of constraints (equality and inequality) +#' the argument is a function cFunc that takes a parameter and returns the +#' evaluated list of constraints, the sparse Jacobian and the constraints type +#' the second argument is an intial parameter value +#' cFunc should return the coloring if called with coloring=TRUE +#' @export +#' @example examples/mpeccable.csolve.ex.r +mpeccable.csolve <- function( cFunc, x0, vars ) { + + # Compute Jacobian Structure + # -------------------------- + + # @todo deal with inequality constraints + + # get the Jacobian from constraints + theta = mpec.vars.collate(x0,vars,coloring=TRUE) + res = cFunc(theta) + + # create the error for each equality constraint, not for inequality + C.MSE = res[['C.MSE']] + + # I create the error and append it + optim.err = FDiff(rep(0,dim(C.MSE)[1]),'optim.err') + C.MSE = C.MSE + optim.err + constraint_lb = 0 + constraint_ub = 0 + + # append optim.err to list of vars + vars[['optim.err']] = dim(optim.err)[1] + + # exctract the sparse structure of the jacobian + eval_jac_g_structure <- make.sparse( C.MSE@J!=0) # needs to + + # Create the objective function and the constraint function for ipopt + # ------------------------------------------------------------------- + + # ------ objective function ------ + eval_f <- function(x,private) { + # extract the squared errors and sum them + ps = mpec.vars.collate(x,private$vars) + return(sum( (ps[['optim.err']]@F)^2)) + } + # ------ objective function gradient ------ + eval_grad_f <- function(x,private) { + # extract the errors and return them + ps = mpec.vars.collate(x,private$vars) + ranges = computeVarRanges(private$vars) + r = mpec.vars.init(private$vars,0) + r[ranges[['optim.rr']]] = ps[['optim.err']]@F + return(r) + } + # ------ constraint function gradient ------ + eval_g <- function( x ,private) { + ps = mpec.vars.collate(x,private$vars) # extract the parameters + optim.err = ps[['optim.err']] # extract the errors + res = private$cFunc(ps) # evaluate the constraints + C.MSE = res[['C.MSE']] + optim.err # add the errors + return(C.MSE@F) + } + # ------ constraint function gradient ------ + eval_jac_g <- function( x ,private) { + ps = mpec.vars.collate(x,private$vars) # extract the parameters + optim.err = ps[['optim.err']] # extract the errors + res = private$cFunc(ps) # evaluate the constraints + C.MSE = res[['C.MSE']] + optim.err # add the errors + + # extract sparse structure of the JAC + res = getSparseValues(C.MSE@J,private$eval_jac_g_structure) + return(res) + } + + x0.augmented = c(x0,optim.err@F) + lb = rep(-Inf,length(x0.augmented)) + ub = rep(Inf,length(x0.augmented)) + + private = list(vars=vars,eval_jac_g_structure=eval_jac_g_structure,cFunc=cFunc) + + length(eval_jac_g(x0.augmented,private)) + length(unlist(eval_jac_g_structure)) + + # Call the optimizer + # ------------------ + res = ipoptr( + x0=x0.augmented, + lb=lb, + ub=ub, + eval_f=eval_f, + eval_grad_f=eval_grad_f, + eval_g=eval_g, + eval_jac_g=eval_jac_g, + eval_jac_g_structure=eval_jac_g_structure, + private=private) + + return(res) +} + + + + + + + diff --git a/R/utils.R b/R/utils.R index 000a42f..f7c4617 100644 --- a/R/utils.R +++ b/R/utils.R @@ -65,4 +65,18 @@ applyColoring <- function(F) { return(F) } +#' variables for the Jacobian columns +#' @export +getSparseValues <- function(M,S) +{ + r = rep(0,length(unlist(S))) + i = 1 + for (row in 1:length(S)) { + for (col in S[[row]]) { + r[i] = M[row,col]; + } + } + + return(r) +} diff --git a/examples/mpeccable.csolve.ex.r b/examples/mpeccable.csolve.ex.r index 93de7ba..4e1ee3c 100644 --- a/examples/mpeccable.csolve.ex.r +++ b/examples/mpeccable.csolve.ex.r @@ -1,5 +1,7 @@ -# The user writes the following: -# ------------------------------ +# This is very simple example that fits 3 splines +# on a very small grid. It is supposed to test +# the piping + require(mpeccable) require(ggplot2) require(ipoptr) @@ -43,8 +45,4 @@ cFunc <- function(params) { return(list(C.MSE=R)) } -# small test -ps = mpec.vars.collate(x0,vars,coloring=TRUE) -res = cFunc(ps) - -res = mpeccable.solve( cFunc, x0, vars ) +res = mpeccable.csolve( cFunc, x0, vars ) diff --git a/test/solve.mpec.r b/test/solve.mpec.r deleted file mode 100644 index 03f5230..0000000 --- a/test/solve.mpec.r +++ /dev/null @@ -1,175 +0,0 @@ - - -# from user point of view -#' solves a set of constraints (equality and inequality) -#' the argument is a function cFunc that takes a parameter and returns the -#' evaluated list of constraints, the sparse Jacobian and the constraints type -#' the second argument is an intial parameter value -#' cFunc should return the coloring if called with coloring=TRUE -#' @export -#' @example examples/mpeccable.csolve.r -mpeccable.csolve <- function( cFunc, x0, vars ) { - - # Compute Jacobian Structure - # -------------------------- - - # @todo bounds on parameters - # @todo deal with inequality constraints - - # get the Jacobian from constraints - theta = mpec.vars.collate(x0,vars,coloring=TRUE) - res = cFunc(ps) - - # create the error for each equality constraint, not for inequality - C.MSE = res[['C.MSE']] - - # I create the error and append it - optim.err = FDiff(rep(0,dim(C.MSE)[1]),'optim.err') - C.MSE = C.MSE + optim.err - constraint_lb = 0 - constraint_ub = 0 - - # append optim.err to list of vars - vars[['optim.err']] = dim(optim.err)[1] - - # exctract the sparse structure of the jacobian - eval_jac_g_structure <- make.sparse( C.MSE@J!=0) # needs to - - # Create the objective function and the constraint function for ipopt - # ------------------------------------------------------------------- - - # ------ objective function ------ - eval_f <- function(x) { - # extract the squared errors and sum them - ps = mpec.vars.collate(x,vars) - return(sum( (ps[['optim.err']]@F)^2)) - } - # ------ objective function gradient ------ - eval_grad_f <- function(x) { - # extract the errors and return them - ps = mpec.vars.collate(x,vars) - return(2*ps[['optim.err']]@F) - } - # ------ constraint function gradient ------ - eval_g <- function( x ) { - ps = mpec.vars.collate(x,vars) # extract the parameters - optim.err = ps[['optim.err']] # extract the errors - res = cFunc(ps) # evaluate the constraints - C.MSE = res[['C.MSE']] + optim.err # add the errors - return(C.MSE@F) - } - # ------ constraint function gradient ------ - eval_jac_g <- function( x ) { - ps = mpec.vars.collate(x,vars) # extract the parameters - optim.err = ps[['optim.err']] # extract the errors - res = cFunc(ps) # evaluate the constraints - C.MSE = res[['C.MSE']] + optim.err # add the errors - - # extract sparse structure of the JAC - return(t(C.MSE@J)@x) - } - - x0.augmented = c(x0,optim.err@F) - - # Call the optimizer - # ------------------ - res = ipoptr( - x0=x0, - eval_f=eval_f, - eval_grad_f=eval_grad_f, - eval_g=eval_g, - eval_jac_g=eval_jac_g, - eval_jac_g_structure=eval_jac_g_structure, - constraint_lb=constraint_lb, - constraint_ub=constraint_ub, - ipoptr_environment = environment()) - - return(res) -} - - - - -# develop a class SolveMpec that wraps the whole thing and sends it off to a solver - -setClass(Class="SolveMpec", representation(Func = "FDiff", x0="numeric", ub="numeric", lb="numeric",sol="list")) - -setMethod("print","SolveMpec", - function(x,...){ - cat("*** Class SolveMpec, print method *** \n") - cat("* initial value = "); print(x@x0) - cat("* ub = "); print(x@ub) - cat("* lb = "); print(x@lb) - cat("* Func = "); print(x@Func) - cat("* Solution = "); print(x@sol) - cat("*** End print (SolveMpec) *** \n") - } - ) - -setMethod(f="show",signature="SolveMpec", - definition=function(object){ - cat("*** Class SolveMpec, show method *** \n") - cat("\n") - cat("* optimization parameters\n") - cat("* -----------------------\n") - cat("* initial value = "); print(head(object@x0,20)) - cat("* ub = "); print(head(object@ub,20)) - cat("* lb = "); print(head(object@lb,20)) - cat("\n") - cat("* Elements of Functional Representation:\n") - cat("* --------------------------------------\n") - cat("* function Values F = \n"); print(head(object@Func@F,20)) - cat("* choice variables = \n"); - print(object@Func@vars) - if (length(object@Func@J)>0){ - cat("* Jacobian[max 30 by 30] = "); - rows <- nrow(object@Func@J) - cols <- ncol(object@Func@J) - print(object@Func@J[1:min(min(rows,cols),30),1:min(min(rows,cols),30)]) - } else { - cat("* empty jacobian slot\n") - } - cat("*** End show (SolveMpec) *** \n") - } - ) - - -# try print and show -M <- new("SolveMpec",Func=R,x0=rnorm(100),ub=rep(3,10),lb=runif(19)) -M - -# initialize an empty object and show it -M2 <- new("SolveMpec") -M2 - - - -# the core function: solve() -setMethod(f="solve",signature="SolveMpec", - definition=function(obj){ - ... - res <- ipoptr(x0=obj@x0, - lb=obj@lb, - ub=obj@ub, - eval_f=objective, - eval_grad_f=grad, - eval_g=constfun, - eval_jac_g=jacfun, - eval_jac_g=constfun, - eval_jac_g_structure=jac.struc) - # you get jac.struc outside of here by - # jac.struc <- make.sparse( R@J ) and this will remain fixed throughout the process. - - # process result - - # store results in respective slots - - # return solver status code - } - ) - - - - - -