ipoptr: R interface to Ipopt

# Copyright (C) 2010 Jelmer Ypma. All Rights Reserved.
# This code is published under the Eclipse Public License.
#
# File:   ipoptr.R
# Author: Jelmer Ypma
# Date:   18 April 2010
#
# Changelog:
#   09/03/2012: Added outputs, z_L, z_U, constraints, lambda (thanks to Michael Schedl)
#   09/03/2012: Removed ipoptr_environment because this caused a bug in combination with 
#               data.table and it wasn't useful (thanks to Florian Oswald for reporting)
#
# Input: 
#		x0 : vector with initial values
#		eval_f : function to evaluate objective function
#		eval_grad_f : function to evaluate gradient of objective function
#		lb : lower bounds of the control
#		ub : upper bounds of the control
#		eval_g : function to evaluate (non-)linear constraints that should hold in the solution
#		eval_jac_g : function to evaluate the jacobian of the (non-)linear constraints that should hold in the solution
#		eval_jac_g_structure : sparseness structure of the jacobian
#		constraint_lb : lower bounds of the (non-)linear constraints
#		constraint_ub : upper bounds of the (non-)linear constraints
#		eval_h : function to evaluate the hessian
#		eval_h_structure : sparseness structure of the hessian
#		opts : list with options that are passed to Ipopt
#       ... : arguments that will be passed to user-defined functions
#
# Output: structure with inputs and
#		call        : the call that was made to solve
#		status      : integer value with the status of the optimization (0 is success)
#		message     : more informative message with the status of the optimization
#		iterations  : number of iterations that were executed
#		objective   : final value of the objective function
#		solution    : final values for the controls
#       z_L         : final values for the lower bound multipliers
#       z_U         : final values for the upper bound multipliers
#       constraints : final values for the constraints
#       lambda      : final values for the Lagrange mulipliers

ipoptr <-
function( x0, 
          eval_f, 
          eval_grad_f, 
          lb = NULL, 
          ub = NULL, 
          eval_g = function( x ) { return( numeric(0) ) }, 
          eval_jac_g = function( x ) { return( numeric(0) ) }, 
          eval_jac_g_structure = list(),
          constraint_lb = numeric(0), 
          constraint_ub = numeric(0),
          eval_h = NULL,
          eval_h_structure = NULL,
          opts = list(),
          ... ) {
    
    # define 'infinite' lower and upper bounds of the control if they haven't been set
    if ( is.null( lb ) ) { lb <- rep( -Inf, length(x0) ) }
    if ( is.null( ub ) ) { ub <- rep(  Inf, length(x0) ) }

    # internal function to check the arguments of the functions
    checkFunctionArguments <- function( fun, arglist, funname ) {
		if( !is.function(fun) ) stop(paste(funname, " must be a function\n", sep = ""))
		
        # determine function arguments
        fargs <- formals(fun)
        
		if ( length(fargs) > 1 ) {
            # determine argument names user-defined function
			argnames_udf <- names(fargs)[2:length(fargs)]	# remove first argument, which is x
            
            # determine argument names that where supplied to ipoptr()
			argnames_supplied <- names(arglist)
            
            # determine which arguments where required but not supplied
			m1 = match(argnames_udf, argnames_supplied)
			if( any(is.na(m1)) ){
				mx1 = which( is.na(m1) )
				for( i in 1:length(mx1) ){
					stop(paste(funname, " requires argument '", argnames_udf[mx1], "' but this has not been passed to the 'ipoptr' function.\n", sep = ""))
				}
			}
            
            # determine which arguments where supplied but not required
			m2 = match(argnames_supplied, argnames_udf)
			if( any(is.na(m2)) ){
				mx2 = which( is.na(m2) )
				for( i in 1:length(mx2) ){
					stop(paste("'", argnames_supplied[mx2], "' passed to (...) in 'ipoptr' but this is not required in the ", funname, " function.\n", sep = ""))
				}
			}
		}
		return( 0 )
	}
    
    # extract list of additional arguments and check user-defined functions
    arglist <- list(...)
	checkFunctionArguments( eval_f, arglist, 'eval_f' )
	checkFunctionArguments( eval_grad_f, arglist, 'eval_grad_f' )
    
    num.constraints <- length( constraint_lb )
    if ( num.constraints > 0 ) {
        checkFunctionArguments( eval_g, arglist, 'eval_g' )
        checkFunctionArguments( eval_jac_g, arglist, 'eval_jac_g' )
    }
    
    # write wrappers around user-defined functions to pass additional arguments
    eval_f_wrapper = function(x){ eval_f(x, ...) }
    eval_grad_f_wrapper = function(x){ eval_grad_f(x, ...) }
    
    if ( num.constraints > 0 ) {
        eval_g_wrapper = function( x ) { eval_g(x, ...) }
        eval_jac_g_wrapper = function( x ) { eval_jac_g(x, ...) }
    } else {
        eval_g_wrapper = function( x ) { return( numeric(0) ) } 
        eval_jac_g_wrapper = function( x ) { return( numeric(0) ) } 
    }

    # approximate Hessian
    if ( is.null( eval_h ) && is.null( eval_h_structure ) ) {
        opts$hessian_approximation <- "limited-memory"
        eval_h_wrapper = NULL
    } else {
        checkFunctionArguments( eval_h, c( arglist, obj_factor=0, hessian_lambda=0 ), 'eval_h' )
        eval_h_wrapper = function( x, obj_factor, hessian_lambda ) { eval_h(x, obj_factor, hessian_lambda, ...) }
    }
	


    # build ipoptr object
    ret <- list( "x0"=x0, 
                 "eval_f"=eval_f_wrapper, 
                 "eval_grad_f"=eval_grad_f_wrapper, 
                 "lower_bounds"=lb, 
                 "upper_bounds"=ub, 
                 "eval_g"=eval_g_wrapper, 
                 "eval_jac_g"=eval_jac_g_wrapper, 
                 "constraint_lower_bounds"=constraint_lb, 
                 "constraint_upper_bounds"=constraint_ub, 
                 "eval_jac_g_structure"=eval_jac_g_structure,
                 "eval_h"=eval_h_wrapper,
                 "eval_h_structure"=eval_h_structure,
                 "options"=get.option.types(opts) )
    
    attr(ret, "class") <- "ipoptr"
    
    # add the current call to the list
    ret$call <- match.call()
    
    # check whether we have a correctly formed ipoptr object
    is.ipoptr( ret )
    
    # pass ipoptr object to C code
    solution <- .Call( IpoptRSolve, ret )
    
    # add solution variables to object
    ret$status      <- solution$status
    ret$message     <- solution$message
    ret$iterations  <- solution$iterations
    ret$objective   <- solution$objective
    ret$solution    <- solution$solution
    ret$z_L         <- solution$z_L
    ret$z_U         <- solution$z_U
    ret$constraints <- solution$constraints
    ret$lambda      <- solution$lambda
    
    return( ret )
}