[med-svn] [r-cran-rsolnp] 06/12: New upstream version 1.16+dfsg
Andreas Tille
tille at debian.org
Wed Nov 29 20:25:49 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-rsolnp.
commit a042818ff83506d068eff6b735e3471b938567ac
Author: Andreas Tille <tille at debian.org>
Date: Wed Nov 29 21:06:05 2017 +0100
New upstream version 1.16+dfsg
---
ChangeLog | 122 +++++++++++
DESCRIPTION | 19 ++
MD5 | 18 ++
NAMESPACE | 8 +
R/benchmarks.R | 492 +++++++++++++++++++++++++++++++++++++++++++
R/gosolnp.R | 459 ++++++++++++++++++++++++++++++++++++++++
R/helpers.R | 351 +++++++++++++++++++++++++++++++
R/solnp.R | 349 +++++++++++++++++++++++++++++++
R/subnp.R | 563 ++++++++++++++++++++++++++++++++++++++++++++++++++
debian/changelog | 14 --
debian/compat | 1 -
debian/control | 24 ---
debian/copyright | 31 ---
debian/rules | 4 -
debian/source/format | 1 -
debian/watch | 3 -
inst/CITATION | 32 +++
inst/COPYRIGHTS | 10 +
inst/doc/manual.pdf | Bin 0 -> 132215 bytes
man/Rsolnp-package.Rd | 25 +++
man/benchmark.Rd | 60 ++++++
man/benchmarkids.Rd | 29 +++
man/gosolnp.Rd | 247 ++++++++++++++++++++++
man/solnp.Rd | 131 ++++++++++++
man/startpars.Rd | 172 +++++++++++++++
25 files changed, 3087 insertions(+), 78 deletions(-)
diff --git a/ChangeLog b/ChangeLog
new file mode 100644
index 0000000..0c213fb
--- /dev/null
+++ b/ChangeLog
@@ -0,0 +1,122 @@
+2015-07-02 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.16
+ * Fix to gosolnp for one parameter case thanks to
+David Lawrence Miller.
+ * Fix to parallel gosolnp not passing multiple function
+arguments (bug reported by Alex Karagiannis).
+ * rep argument in solnp.R had length rather than length.out (leading
+ to cases name conflict).
+ * Fix to pass new CRAN checks which only attach base.
+
+2013-04-10 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.15
+ * Fixes to pass 3.0.0
+
+2012-12-05 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.14
+ * Change to parallel functionality for gosolnp. Now makes use of the
+ parallel package where the user has to pass an initialized cluster
+ object (also fixes Windows paralllel functionality problems).
+
+2012-05-22 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.13
+ * Fix to rtruncnorm call.
+
+2012-05-22 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.12
+ * Fix to .onLoad in zzz.R
+
+2011-07-15 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.11
+ * Added a penalty barrier function for use in the gosolnp random sampling and evaluation.
+ * New function startpars directly returns a set of starting parameters ranked according
+ to either a direct evaluation of the objective (excluding any inequality violations) else
+ a penalty barrier function comprising the objective and constraints.
+
+2011-07-06 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.1
+ * Correction to bug when terminating without solution (introduced in last update).
+
+2011-05-19 Alexios Ghalanos <alexios at 4dscape.com>
+ * DESCRIPTION (Version): New version is 1.0-9
+ * Added more traps to catch errors in inverting the hessian during the inner step.
+ Convergence codes are now (0 - solution), (1 - maximum iterations without tolerance),
+ (2 - no convergence, failure to invert Hessian).
+
+2011-01-03 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-8
+ * Parallel functionality changed in 'gosolnp' function to make use of both multicore and snowfall
+ packages (i.e. windows multicore functionality now possible). This also opens the possibility for
+ future extensions to use snowfall with Rmpi for cluster processing.
+
+2010-10-01 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-7
+ * Some additional clarification in the documentation on the control parameters.
+ * Corrections to gosolnp checks on the inequality violations during random parameter generation
+ in the case of one-dimensional constraints (thanks to Matthieu Stigler).
+ * Convergence checks in the gosolnp with multiple restarts method to avoid exiting on failure from
+ solnp.
+
+2010-09-22 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-6
+ * Correction to gosolnp bugs introduced in 1.0-5.
+ * Parameter vector (pars) now retains names across calls to fun, eqfun and ineqfun
+ (previously these were removed).
+
+2010-09-03 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-5
+ * Trace parameter was not correctly passed to gosolnp. Also minor
+ corrections to trace parameter in solnp (i.e. warnings off).
+
+2010-06-27 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-4
+ * Corrections to subnp function when single parameter passed
+ to solver (effectively added dimension arguments to 'diag'
+ function when updating the Hessian).
+
+2010-03-20 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0-3
+ * Corrections to documentation
+
+2010-03-10 Stefan Theussl <stefan.theussl at wu.ac.at>
+
+ * DESCRIPTION (Version): New version is 1.0-2
+ * Correction to default solver control parameters
+
+
+2010-03-05 Stefan Theussl <stefan.theussl at wu.ac.at>
+
+ * DESCRIPTION (Version): New version is 1.0-1
+ * CITATION file: added
+
+2010-03-10 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 1.0
+ * Release version for CRAN
+ * Added a function to generate random starting parameters and allow
+ for multiple restarts of the solver.
+ * Some more benchmark problems to test the gosolnp function.
+
+2009-09-15 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 0.3.
+ * Changes to function inputs and benchmark suite and method
+ of function evaluation.
+ * Minor Fixes
+ * More Checks to input functions
+
+2009-05-21 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 0.2.
+ * Benchmark Problem suite created with more problems added.
+
+2009-05-14 Alexios Ghalanos <alexios at 4dscape.com>
+
+ * DESCRIPTION (Version): New version is 0.1.
+ * First upload of Rsolnp.
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..7ab1e4a
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,19 @@
+Package: Rsolnp
+Type: Package
+Title: General Non-Linear Optimization
+Version: 1.16
+Date: 2015-07-02
+Author: Alexios Ghalanos and Stefan Theussl
+Maintainer: Alexios Ghalanos <alexios at 4dscape.com>
+Depends: R (>= 2.10.0)
+Imports: truncnorm, parallel, stats
+Description: General Non-linear Optimization Using Augmented Lagrange Multiplier Method.
+LazyLoad: yes
+License: GPL
+Repository: CRAN
+Repository/R-Forge/Project: rino
+Repository/R-Forge/Revision: 102
+Repository/R-Forge/DateTimeStamp: 2015-07-04 02:08:32
+Date/Publication: 2015-12-28 09:01:11
+NeedsCompilation: no
+Packaged: 2015-07-08 12:20:39 UTC; rforge
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..a720b61
--- /dev/null
+++ b/MD5
@@ -0,0 +1,18 @@
+ff7e87808b9231c0f9c9c9f6e3de8bd0 *ChangeLog
+1b9fd85920581411d63c4ac5e4bb1fc5 *DESCRIPTION
+fe2626ab83911737139dde60a18e2c47 *NAMESPACE
+a6419b2b559df95a04e004d6f2af4e1c *R/benchmarks.R
+995f7c1149d4fd5663c3827c21e4209c *R/gosolnp.R
+9a31d8d6bddac7e89cc65861bc8e43b2 *R/helpers.R
+b34a820c2cbbb1879eb6eff9b5eeb3dd *R/solnp.R
+bb03327b209d67b2a2885eebc6394528 *R/subnp.R
+b21f2dbd4bb3f73936ebd790600bbc81 *inst/CITATION
+8fd397e50691113d3d9befbd6964654f *inst/COPYRIGHTS
+bc76d8b2815126d9d9e606a1fafd3e63 *inst/doc/manual.pdf
+31458a798fe8f9ff6bc6ee9731ee0e62 *inst/doc/manual.ps
+e33c76896f5b19bfca547e82bcfc5454 *man/Rsolnp-package.Rd
+b0a3135e881d767e2e38b36facccb85e *man/benchmark.Rd
+188935cba607ac4833ddde85e2bf59c5 *man/benchmarkids.Rd
+7e1a6fcbba6cfe7f32df4b548888564b *man/gosolnp.Rd
+e577612e27bec990ec5a5a90159b84bc *man/solnp.Rd
+4d4f57a436e792b0a136b8b03a47ff39 *man/startpars.Rd
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..b94459a
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,8 @@
+import(truncnorm)
+importFrom("stats", "rnorm", "runif")
+importFrom(parallel, parLapply, clusterExport, clusterEvalQ)
+export("solnp")
+export("gosolnp")
+export("startpars")
+export("benchmark")
+export("benchmarkids")
diff --git a/R/benchmarks.R b/R/benchmarks.R
new file mode 100644
index 0000000..517a607
--- /dev/null
+++ b/R/benchmarks.R
@@ -0,0 +1,492 @@
+#################################################################################
+##
+## R package Rsolnp by Alexios Ghalanos and Stefan Theussl Copyright (C) 2009-2013
+## This file is part of the R package Rsolnp.
+##
+## The R package Rsolnp is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## The R package Rsolnp is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+#################################################################################
+
+benchmarkids <- function()
+{
+return(c("Powell", "Wright4", "Wright9", "Alkylation", "Entropy", "Box", "RosenSuzuki",
+ "Electron", "Permutation"))
+}
+
+
+benchmark <- function( id = "Powell")
+{
+ if( !any(benchmarkids() == id[ 1L ]) )
+ stop( "invalid benchmark id" )
+ ans = switch(id,
+ Powell = .powell(),
+ Wright4 = .wright4(),
+ Wright9 = .wright9(),
+ Alkylation = .alkylation(),
+ Entropy = .entropy(),
+ Box = .box(),
+ RosenSuzuki = .rosensuzuki(),
+ Electron = .electron(),
+ Permutation = .permutation())
+ return(ans)
+}
+
+.powell = function()
+{
+ .fn1 = function(x)
+ {
+ exp(x[1]*x[2]*x[3]*x[4]*x[5])
+ }
+
+ .eqn1 = function(x){
+ z1=x[1]*x[1]+x[2]*x[2]+x[3]*x[3]+x[4]*x[4]+x[5]*x[5]
+ z2=x[2]*x[3]-5*x[4]*x[5]
+ z3=x[1]*x[1]*x[1]+x[2]*x[2]*x[2]
+ return(c(z1,z2,z3))
+ }
+
+ .x0 = c(-2, 2, 2, -1, -1)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, eqfun = .eqn1, eqB = c(10,0,-1), control=ctrl)
+ minos = list()
+ minos$fn = 0.05394985
+ minos$pars = c(-1.717144, 1.595710, 1.827245, 0.763643, 0.763643)
+ minos$nfun = 524
+ minos$iter = 12
+ minos$elapsed = 0.2184
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("Powell's exponential problem is a function of five variables with
+three nonlinear equality constraints on the variables.")
+
+ return(bt)
+}
+
+.wright4 = function()
+{
+ .fn1 = function(x)
+ {
+ (x[1]-1)^2+(x[1]-x[2])^2+(x[2]-x[3])^3+(x[3]-x[4])^4+(x[4]-x[5])^4
+ }
+
+ .eqn1 = function(x){
+ z1=x[1]+x[2]*x[2]+x[3]*x[3]*x[3]
+ z2=x[2]-x[3]*x[3]+x[4]
+ z3=x[1]*x[5]
+ return(c(z1,z2,z3))
+ }
+
+ .x0 = c(1, 1, 1, 1, 1)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, eqfun = .eqn1, eqB = c(2+3*sqrt(2),-2+2*sqrt(2),2), control=ctrl)
+ minos = list()
+ minos$fn = 0.02931083
+ minos$pars = c(1.116635, 1.220442, 1.537785, 1.972769, 1.791096)
+ minos$nfun = 560
+ minos$iter = 9
+ minos$elapsed = 0.249
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("Wright's fourth problem is a function of five variables with
+three non linear equality constraints on the variables. This popular
+test problem has several local solutions and taken from Wright (1976).")
+ return(bt)
+}
+
+.wright9 = function()
+{
+ .fn1 = function(x)
+ {
+ 10*x[1]*x[4]-6*x[3]*x[2]*x[2]+x[2]*(x[1]*x[1]*x[1])+
+ 9*sin(x[5]-x[3])+x[5]^4*x[4]*x[4]*x[2]*x[2]*x[2]
+ }
+
+ .ineqn1 = function(x){
+ z1=x[1]*x[1]+x[2]*x[2]+x[3]*x[3]+x[4]*x[4]+x[5]*x[5]
+ z2=x[1]*x[1]*x[3]-x[4]*x[5]
+ z3=x[2]*x[2]*x[4]+10*x[1]*x[5]
+ return(c(z1,z2,z3))
+ }
+ ineqLB = c(-100, -2, 5)
+ ineqUB = c(20, 100, 100)
+ .x0 = c(1, 1, 1, 1, 1)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, ineqfun = .ineqn1, ineqLB = ineqLB, ineqUB = ineqUB, control=ctrl)
+ minos = list()
+ minos$fn = -210.4078
+ minos$pars = c(-0.08145219, 3.69237756, 2.48741102, 0.37713392, 0.17398257)
+ minos$nfun = 794
+ minos$iter = 11
+ minos$elapsed = 0.281
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("Wright's ninth problem is a function of five variables with three
+non linear inequality constraints on the variables. This popular test
+problem has several local solutions and taken from Wright (1976).")
+ return(bt)
+}
+
+.alkylation = function()
+{
+ .fn1 = function(x)
+ {
+ -0.63*x[4]*x[7]+50.4*x[1]+3.5*x[2]+x[3]+33.6*x[5]
+ }
+
+ .eqn1 = function(x){
+ z1=98*x[3]-0.1*x[4]*x[6]*x[9]-x[3]*x[6]
+ z2=1000*x[2]+100*x[5]-100*x[1]*x[8]
+ z3=122*x[4]-100*x[1]-100*x[5]
+ return(c(z1,z2,z3))
+ }
+ .ineqn1 = function(x){
+ z1=(1.12*x[1]+0.13167*x[1]*x[8]-0.00667*x[1]*x[8]*x[8])/x[4]
+ z2=(1.098*x[8]-0.038*x[8]*x[8]+0.325*x[6]+57.25)/x[7]
+ z3=(-0.222*x[10]+35.82)/x[9]
+ z4=(3*x[7]-133)/x[10]
+ return(c(z1,z2,z3,z4))
+ }
+ ineqLB = c(0.99,0.99,0.9,0.99)
+ ineqUB = c(100/99,100/99,10/9,100/99)
+ eqB = c(0,0,0)
+ LB = c(0,0,0,10,0,85,10,3,1,145)
+ UB = c(20,16,120,50,20,93,95,12,4,162)
+ .x0 = c(17.45,12,110,30,19.74,89.2,92.8,8,3.6,155)
+ ctrl = list(rho = 0, trace=0)
+ ans = solnp(.x0, fun = .fn1, eqfun = .eqn1, eqB = eqB, ineqfun = .ineqn1, ineqLB = ineqLB,
+ ineqUB = ineqUB, LB = LB, UB = UB, control = ctrl)
+ minos = list()
+ minos$fn = -172.642
+ minos$pars = c(16.996427, 16.000000, 57.685751, 30.324940, 20.000000, 90.565147, 95.000000, 10.590461, 1.561636, 153.535354)
+ minos$nfun = 2587
+ minos$iter = 13
+ minos$elapsed = 0.811
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("The Alkylation problem models a simplified alkylation process. It
+is a function of ten variables with four non linear inequality and
+three non linear equality constraints as well as variable bounds.
+The problem is taken from Locke and Westerberg (1980).")
+ return(bt)
+}
+
+.entropy = function()
+{
+ .fn1 = function(x)
+ {
+ m = length(x)
+ f = 0
+ for(i in 1:m){
+ f = f-log(x[i])
+ }
+ ans = f-log(.vnorm(x-1) + 0.1)
+ ans
+ }
+
+ .eqn1 = function(x){
+ sum(x)
+ }
+ eqB = 10
+ LB = rep(0,10)
+ UB = rep(1000,10)
+ .x0 = runif(10, 0, 1000)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, eqfun = .eqn1, eqB = eqB, LB = LB, UB = UB, control=ctrl)
+ minos = list()
+ minos$fn = 0.1854782
+ minos$pars = c(2.2801555, 0.8577605, 0.8577605, 0.8577605, 0.8577605, 0.8577605,
+ 0.8577605, 0.8577605, 0.8577605, 0.8577605)
+ minos$nfun = 886
+ minos$iter = 4
+ minos$elapsed = 0.296
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("The Entropy problem is non convex in n variables with one linear
+equality constraint and variable positivity bounds.")
+ return(bt)
+}
+
+.box = function()
+{
+ .fn1 = function(x)
+ {
+ -x[1]*x[2]*x[3]
+ }
+
+ .eqn1 = function(x){
+ 4*x[1]*x[2]+2*x[2]*x[3]+2*x[3]*x[1]
+ }
+
+ eqB = 100
+ LB = rep(1, 3)
+ UB = rep(10, 3)
+
+ .x0 = c(1.1, 1.1, 9)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, eqfun = .eqn1, eqB = eqB, LB = LB, UB = UB, control=ctrl)
+ minos = list()
+ minos$fn = -48.11252
+ minos$pars = c(2.886751, 2.886751, 5.773503)
+ minos$nfun = 394
+ minos$iter = 9
+ minos$elapsed = 0.156
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("The box problem is a function of three variables with one non
+linear equality constraint and variable bounds.")
+ return(bt)
+}
+
+.rosensuzuki = function()
+{
+ .fn1 = function(x)
+ {
+ x[1]*x[1]+x[2]*x[2]+2*x[3]*x[3]+x[4]*x[4]-5*x[1]-5*x[2]-21*x[3]+7*x[4]
+ }
+
+ .ineqn1 = function(x){
+ z1=8-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[4]*x[4]-x[1]+x[2]-x[3]+x[4]
+ z2=10-x[1]*x[1]-2*x[2]*x[2]-x[3]*x[3]-2*x[4]*x[4]+x[1]+x[4]
+ z3=5-2*x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-2*x[1]+x[2]+x[4]
+ return(c(z1,z2,z3))
+ }
+ ineqLB = rep(0, 3)
+ ineqUB = rep(1000, 3)
+ .x0 = c(1, 1, 1, 1)
+ ctrl=list(trace=0)
+ ans = solnp(.x0, fun = .fn1, ineqfun = .ineqn1, ineqLB = ineqLB, ineqUB = ineqUB, control=ctrl)
+ minos = list()
+ minos$fn = -44
+ minos$pars = c(2.502771e-07, 9.999997e-01, 2.000000e+00, -1.000000e+00)
+ minos$nfun = 527
+ minos$iter = 12
+ minos$elapsed = 0.203
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ minos = rbind(round(minos$fn, 5L),
+ round(minos$iter, 0L),
+ round(0, 0L),
+ round(minos$nfun, 0L),
+ round(minos$elapsed, 3L),
+ matrix(round(minos$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ attr(bt, "description") = paste("The Rosen-Suzuki problem is a function of four variables with
+three nonlinear inequality constraints on the variables. It is taken
+from Problem 43 of Hock and Schittkowski (1981).")
+ return(bt)
+}
+
+
+
+#----------------------------------------------------------------------------------
+# Some Problems in Global Optimization
+#----------------------------------------------------------------------------------
+
+
+# Distribution of Electrons on a Sphere
+# Given n electrons, find the equilibrium state distribution (of minimal Coulomb potential)
+# of the electrons positioned on a conducting sphere. This model is from the COPS benchmarking suite.
+# See http://www-unix.mcs.anl.gov/~more/cops/.
+
+.electron = function()
+{
+ gofn = function(dat, n)
+ {
+
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
+ jj = matrix(1:n, ncol = n, nrow = n)
+ ij = which(ii<jj, arr.ind = TRUE)
+ i = ij[,1]
+ j = ij[,2]
+ # Coulomb potential
+ potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
+ potential
+ }
+
+ goeqfn = function(dat, n)
+ {
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ apply(cbind(x^2, y^2, z^2), 1, "sum")
+ }
+
+ n = 25
+ LB = rep(-1, 3*n)
+ UB = rep(1, 3*n)
+ eqB = rep(1, n)
+ ans = gosolnp(pars = NULL, fixed = NULL, fun = gofn, eqfun = goeqfn, eqB = eqB, LB = LB, UB = UB,
+ control = list(), distr = rep(1, length(LB)), distr.opt = list(outer.iter = 10, trace = 1),
+ n.restarts = 2, n.sim = 20000, rseed = 443, n = 25)
+
+ conopt = list()
+ conopt$fn = 243.813
+ conopt$iter = 33
+ conopt$nfun = NA
+ conopt$elapsed = 0.041
+ conopt$pars = c(-0.0117133872042326, 0.627138691757704, -0.471025867741051, -0.164419761338935, -0.0315460712487934,
+ -0.12718981058582, -0.540049624346613, 0.600346770449059, 0.29796281847713, -0.740960572770077, 0.972512148478245,
+ -0.870858895858346, 0.84178885636396, -0.182471994739506, 0.603293664844919, 0.0834172554171806, 0.51317309921937,
+ 0.260639996237799, -0.0972877803105543, -0.979882559381314, -0.64809471648373, -0.722351411610064, 0.847184430059647,
+ 0.514683899757428, -0.574607272207711, 0.114609211815613, -0.748168886860133, -0.379763612890494, 0.743271243797936,
+ 0.846784034469756, 0.220425955966718, 0.839147591392778, -0.613810163641104, -0.499794531840362, -0.199680552951248,
+ -0.105141937435843, -0.434753057357539, -0.127562956191463, -0.895691740627038, 0.574257349984438, -0.967631920158332,
+ -0.0243647398873149, 0.959445715727407, -0.517406241891138, 0.197677956191858, 0.503867654081605, 0.286619450971711,
+ 0.522509289901788, 0.474911361600545, -0.768915699421274, -0.993341595387613, -0.21665728244143, -0.796187308516752, -0.648470508368975,
+ 0.531000606737778, 0.967075565826839, -0.0646353084836963, 0.512660548728168, -0.813279524362711, 0.641174786133487,
+ 0.207762590603952, -0.229335044471304, -0.524518077390237, -0.405512363446903, 0.55341236880543, 0.23818066376044,
+ 0.857939234263016, -0.107381147942665, 0.850191665834437, 0.0274516931378701, 0.571043453369507, -0.629331175510656,
+ -0.0962423162171412, -0.713834491989006, 0.280348229724225)
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ conopt = rbind(round(conopt$fn, 5L),
+ round(conopt$iter, 0L),
+ round(0, 0L),
+ round(conopt$nfun, 0L),
+ round(conopt$elapsed, 3L),
+ matrix(round(conopt$pars, 5L), ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ colnames(bt) = c("solnp", "conopt")
+ attr(bt, "description") = paste("The equilibrium state distribution (of minimal Coulomb potential)\n of the electrons positioned on a conducting sphere.")
+ return(bt)
+}
+
+# Permutation Problem -- Unique Solution f(x) = 0 and x(i) = i
+
+.permutation = function()
+{
+ .perm = function(x, n, b){
+ F = 0
+ for(k in 1:n){
+ S = 0
+ for(i in 1:n){
+ S = S + ( ( (i^k) + b ) * (( x[i]/i )^k -1))
+ }
+ F = F + S^2
+ }
+ F
+ }
+
+ ans = gosolnp(pars = NULL, fixed = NULL, fun = .perm, eqfun = NULL, eqB = NULL, LB = rep(-4, 4), UB = rep(4, 4),
+ control = list(outer.iter = 25, trace = 1, tol = 1e-9), distr = rep(1, 4), distr.opt = list(),
+ n.restarts = 6, n.sim = 20000, rseed = 99, n = 4, b =0.5)
+
+
+ bt = data.frame( solnp = rbind(round(ans$values[length(ans$values)], 5L),
+ round(ans$outer.iter, 0L),
+ round(ans$convergence, 0L),
+ round(ans$nfuneval, 0L),
+ round(ans$elapsed, 3L),
+ matrix(round(ans$pars, 5L), ncol = 1L)),
+ actual = rbind(0,
+ NA,
+ round(0, 0L),
+ NA,
+ NA,
+ matrix(1:4, ncol = 1L)) )
+ rownames(bt) <- c("funcValue", "majorIter", "exitFlag", "nfunEval", "time(sec)",
+ paste("par.", 1L:length(ans$pars), sep = "") )
+ colnames(bt) = c("solnp", "expected")
+ attr(bt, "description") = paste("Permutation Problem PERM(4,0.5).")
+
+}
\ No newline at end of file
diff --git a/R/gosolnp.R b/R/gosolnp.R
new file mode 100644
index 0000000..a83afd8
--- /dev/null
+++ b/R/gosolnp.R
@@ -0,0 +1,459 @@
+#################################################################################
+##
+## R package Rsolnp by Alexios Ghalanos and Stefan Theussl Copyright (C) 2009-2013
+## This file is part of the R package Rsolnp.
+##
+## The R package Rsolnp is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## The R package Rsolnp is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+#################################################################################
+
+#---------------------------------------------------------------------------------
+# optimization by randomly restarted parameters using simulated parameter strategy
+# Alexios Ghalanos 2010
+#---------------------------------------------------------------------------------
+
+# allowed distributions:
+# 1: uniform (no confidence in the location of the parameter...somewhere in LB-UB space)
+# 2: truncnorm (high confidence in the location of the parameter)
+# 3: normal (Uncertainty in Lower and Upper bounds, but some idea about the dispersion about the location)
+# ...
+
+gosolnp = function(pars = NULL, fixed = NULL, fun, eqfun = NULL, eqB = NULL, ineqfun = NULL, ineqLB = NULL,
+ ineqUB = NULL, LB = NULL, UB = NULL, control = list(), distr = rep(1, length(LB)), distr.opt = list(),
+ n.restarts = 1, n.sim = 20000, cluster = NULL, rseed = NULL, ...)
+{
+ if( !is.null(pars) ) gosolnp_parnames = names(pars) else gosolnp_parnames = NULL
+ if(is.null(control$trace)) trace = FALSE else trace = as.logical(control$trace)
+ if(is.null(control$eval.type)) parmethod = 1 else parmethod = as.integer(min(abs(control$eval.type),2))
+ if(parmethod == 0) parmethod = 1
+ control$eval.type = NULL
+ # use a seed to initialize random no. generation
+ if(is.null(rseed)) rseed = as.numeric(Sys.time()) else rseed = as.integer(rseed)
+ # function requires both upper and lower bounds
+ if(is.null(LB))
+ stop("\ngosolnp-->error: the function requires lower parameter bounds\n", call. = FALSE)
+ if(is.null(UB))
+ stop("\ngosolnp-->error: the function requires upper parameter bounds\n", call. = FALSE)
+ # allow for fixed parameters (i.e. non randomly chosen), but require pars vector in that case
+ if(!is.null(fixed) && is.null(pars))
+ stop("\ngosolnp-->error: you need to provide a pars vector if using the fixed option\n", call. = FALSE)
+ if(!is.null(pars)) n = length(pars) else n = length(LB)
+
+ np = 1:n
+
+ if(!is.null(fixed)){
+ # make unique
+ fixed = unique(fixed)
+ # check for violations in indices
+ if(any(is.na(match(fixed, np))))
+ stop("\ngosolnp-->error: fixed indices out of bounds\n", call. = FALSE)
+ }
+ # check distribution options
+ # truncated normal
+ if(any(distr == 2)){
+ d2 = which(distr == 2)
+ for(i in 1:length(d2)) {
+ if(is.null(distr.opt[[d2[i]]]$mean))
+ stop(paste("\ngosolnp-->error: distr.opt[[,",d2[i],"]] missing mean\n", sep = ""), call. = FALSE)
+ if(is.null(distr.opt[[d2[i]]]$sd))
+ stop(paste("\ngosolnp-->error: distr.opt[[,",d2[i],"]] missing sd\n", sep = ""), call. = FALSE)
+ }
+ }
+ # normal
+ if(any(distr == 3)){
+ d3 = which(distr == 3)
+ for(i in 1:length(d3)) {
+ if(is.null(distr.opt[[d3[i]]]$mean))
+ stop(paste("\ngosolnp-->error: distr.opt[[,",d3[i],"]] missing mean\n", sep = ""), call. = FALSE)
+ if(is.null(distr.opt[[d3[i]]]$sd))
+ stop(paste("\ngosolnp-->error: distr.opt[[,",d3[i],"]] missing sd\n", sep = ""), call. = FALSE)
+ }
+ }
+ # setup cluster exports:
+ if( !is.null(cluster) ){
+ clusterExport(cluster, c("gosolnp_parnames", "fun", "eqfun",
+ "eqB", "ineqfun", "ineqLB", "ineqUB", "LB", "UB"), envir = environment())
+ if( !is.null(names(list(...))) ){
+ # evaluate promises
+ xl = names(list(...))
+ for(i in 1:length(xl)){
+ eval(parse(text=paste(xl[i],"=list(...)[[i]]",sep="")))
+ }
+ clusterExport(cluster, names(list(...)), envir = environment())
+ }
+ clusterEvalQ(cluster, require(Rsolnp))
+ }
+ # initiate random search
+ gosolnp_rndpars = switch(parmethod,
+ .randpars(pars = pars, fixed = fixed, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun, ineqLB = ineqLB,
+ ineqUB = ineqUB, LB = LB, UB = UB, distr = distr,
+ distr.opt = distr.opt, n.restarts = n.restarts,
+ n.sim = n.sim, trace = trace, rseed = rseed,
+ gosolnp_parnames = gosolnp_parnames, cluster = cluster, ...),
+ .randpars2(pars = pars, fixed = fixed, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun, ineqLB = ineqLB,
+ ineqUB = ineqUB, LB = LB, UB = UB, distr = distr,
+ distr.opt = distr.opt, n.restarts = n.restarts,
+ n.sim = n.sim, rseed = rseed, trace = trace,
+ gosolnp_parnames = gosolnp_parnames, cluster = cluster, ...))
+
+ gosolnp_rndpars = gosolnp_rndpars[,1:n, drop = FALSE]
+ # initiate solver restarts
+ if( trace ) cat("\ngosolnp-->Starting Solver\n")
+ solution = vector(mode = "list", length = n.restarts)
+ if( !is.null(cluster) )
+ {
+ clusterExport(cluster, c("gosolnp_rndpars"), envir = environment())
+ solution = parLapply(cluster, as.list(1:n.restarts), fun = function(i) {
+ xx = gosolnp_rndpars[i,]
+ names(xx) = gosolnp_parnames
+ ans = try(solnp(pars = xx, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun,
+ ineqLB = ineqLB, ineqUB = ineqUB,
+ LB = LB, UB = UB,
+ control = control, ...), silent = TRUE)
+ if(inherits(ans, "try-error")){
+ ans = list()
+ ans$values = 1e10
+ ans$convergence = 0
+ ans$pars = rep(NA, length(xx))
+ }
+ return( ans )
+ })
+ } else {
+ solution = lapply(as.list(1:n.restarts), FUN = function(i){
+ xx = gosolnp_rndpars[i,]
+ names(xx) = gosolnp_parnames
+ ans = try(solnp(pars = xx, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun,
+ ineqLB = ineqLB, ineqUB = ineqUB,
+ LB = LB, UB = UB,
+ control = control, ...), silent = TRUE)
+ if(inherits(ans, "try-error")){
+ ans = list()
+ ans$values = 1e10
+ ans$convergence = 0
+ ans$pars = rep(NA, length(xx))
+ }
+ return( ans )
+ })
+ }
+ if(n.restarts>1){
+ best = sapply(solution, FUN = function(x) if(x$convergence!=0) NA else x$values[length(x$values)])
+ if(all(is.na(best)))
+ stop("\ngosolnp-->Could not find a feasible starting point...exiting\n", call. = FALSE)
+ nb = which(best == min(best, na.rm = TRUE))[1]
+ solution = solution[[nb]]
+ if( trace ) cat("\ngosolnp-->Done!\n")
+ solution$start.pars = gosolnp_rndpars[nb,]
+ names(solution$start.pars) = gosolnp_parnames
+ solution$rseed = rseed
+ } else{
+ solution = solution[[1]]
+ solution$start.pars = gosolnp_rndpars[1,]
+ names(solution$start.pars) = gosolnp_parnames
+ solution$rseed = rseed
+ }
+ return(solution)
+}
+
+
+
+startpars = function(pars = NULL, fixed = NULL, fun, eqfun = NULL, eqB = NULL,
+ ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, LB = NULL, UB = NULL,
+ distr = rep(1, length(LB)), distr.opt = list(), n.sim = 20000, cluster = NULL,
+ rseed = NULL, bestN = 15, eval.type = 1, trace = FALSE, ...)
+{
+ if( !is.null(pars) ) gosolnp_parnames = names(pars) else gosolnp_parnames = NULL
+ if(is.null(eval.type)) parmethod = 1 else parmethod = as.integer(min(abs(eval.type),2))
+ if(parmethod == 0) parmethod = 1
+ eval.type = NULL
+ #trace = FALSE
+ # use a seed to initialize random no. generation
+ if(is.null(rseed)) rseed = as.numeric(Sys.time()) else rseed = as.integer(rseed)
+ # function requires both upper and lower bounds
+ if(is.null(LB))
+ stop("\nstartpars-->error: the function requires lower parameter bounds\n", call. = FALSE)
+ if(is.null(UB))
+ stop("\nstartpars-->error: the function requires upper parameter bounds\n", call. = FALSE)
+
+ # allow for fixed parameters (i.e. non randomly chosen), but require pars vector in that case
+ if(!is.null(fixed) && is.null(pars))
+ stop("\nstartpars-->error: you need to provide a pars vector if using the fixed option\n", call. = FALSE)
+ if(!is.null(pars)) n = length(pars) else n = length(LB)
+
+ np = seq_len(n)
+
+ if(!is.null(fixed)){
+ # make unique
+ fixed = unique(fixed)
+ # check for violations in indices
+ if(any(is.na(match(fixed, np))))
+ stop("\nstartpars-->error: fixed indices out of bounds\n", call. = FALSE)
+ }
+
+ # check distribution options
+ # truncated normal
+ if(any(distr == 2)){
+ d2 = which(distr == 2)
+ for(i in 1:length(d2)) {
+ if(is.null(distr.opt[[d2[i]]]$mean))
+ stop(paste("\nstartpars-->error: distr.opt[[,",d2[i],"]] missing mean\n", sep = ""), call. = FALSE)
+ if(is.null(distr.opt[[d2[i]]]$sd))
+ stop(paste("\nstartpars-->error: distr.opt[[,",d2[i],"]] missing sd\n", sep = ""), call. = FALSE)
+ }
+ }
+ # normal
+ if(any(distr == 3)){
+ d3 = which(distr == 3)
+ for(i in 1:length(d3)) {
+ if(is.null(distr.opt[[d3[i]]]$mean))
+ stop(paste("\nstartpars-->error: distr.opt[[,",d3[i],"]] missing mean\n", sep = ""), call. = FALSE)
+ if(is.null(distr.opt[[d3[i]]]$sd))
+ stop(paste("\nstartpars-->error: distr.opt[[,",d3[i],"]] missing sd\n", sep = ""), call. = FALSE)
+ }
+ }
+
+ # setup cluster exports:
+ if( !is.null(cluster) ){
+ clusterExport(cluster, c("gosolnp_parnames", "fun", "eqfun",
+ "eqB", "ineqfun", "ineqLB", "ineqUB", "LB", "UB"), envir = environment())
+ if( !is.null(names(list(...))) ){
+ # evaluate promises
+ xl = names(list(...))
+ for(i in 1:length(xl)){
+ eval(parse(text = paste(xl[i], "=list(...)", "[[" , i, "]]", sep = "")))
+ }
+ clusterExport(cluster, names(list(...)), envir = environment())
+ }
+ if( !is.null(names(list(...))) ) parallel::clusterExport(cluster, names(list(...)), envir = environment())
+ clusterEvalQ(cluster, require(Rsolnp))
+ }
+
+ # initiate random search
+ gosolnp_rndpars = switch(parmethod,
+ .randpars(pars = pars, fixed = fixed, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun, ineqLB = ineqLB, ineqUB = ineqUB,
+ LB = LB, UB = UB, distr = distr, distr.opt = distr.opt,
+ n.restarts = as.integer(bestN), n.sim = n.sim, trace = trace,
+ rseed = rseed, gosolnp_parnames = gosolnp_parnames,
+ cluster = cluster, ...),
+ .randpars2(pars = pars, fixed = fixed, fun = fun, eqfun = eqfun,
+ eqB = eqB, ineqfun = ineqfun, ineqLB = ineqLB, ineqUB = ineqUB,
+ LB = LB, UB = UB, distr = distr, distr.opt = distr.opt,
+ n.restarts = as.integer(bestN), n.sim = n.sim, trace = trace,
+ rseed = rseed, gosolnp_parnames = gosolnp_parnames,
+ cluster = cluster, ...))
+ return(gosolnp_rndpars)
+}
+
+
+.randpars = function(pars, fixed, fun, eqfun, eqB, ineqfun, ineqLB, ineqUB,
+ LB, UB, distr, distr.opt, n.restarts, n.sim, trace = TRUE, rseed,
+ gosolnp_parnames, cluster, ...)
+{
+ if( trace ) cat("\nCalculating Random Initialization Parameters...")
+ N = length(LB)
+ gosolnp_rndpars = matrix(NA, ncol = N, nrow = n.sim * n.restarts)
+ if(!is.null(fixed)) for(i in 1:length(fixed)) gosolnp_rndpars[,fixed[i]] = pars[fixed[i]]
+ nf = 1:N
+ if(!is.null(fixed)) nf = nf[-c(fixed)]
+ m = length(nf)
+ set.seed(rseed)
+ for(i in 1:m){
+ j = nf[i]
+ gosolnp_rndpars[,j] = switch(distr[j],
+ .distr1(LB[j], UB[j], n.restarts*n.sim),
+ .distr2(LB[j], UB[j], n.restarts*n.sim, mean = distr.opt[[j]]$mean, sd = distr.opt[[j]]$sd),
+ .distr3(n.restarts*n.sim, mean = distr.opt[[j]]$mean, sd = distr.opt[[j]]$sd)
+ )
+ }
+
+ if( trace ) cat("ok!\n")
+
+ if(!is.null(ineqfun)){
+ if( trace ) cat("\nExcluding Inequality Violations...\n")
+ ineqv = matrix(NA, ncol = length(ineqLB), nrow = n.restarts*n.sim)
+ # ineqv = t(apply(rndpars, 1, FUN = function(x) ineqfun(x)))
+ if(length(ineqLB) == 1){
+ ineqv = apply(gosolnp_rndpars, 1, FUN = function(x){
+ names(x) = gosolnp_parnames
+ ineqfun(x, ...)} )
+ lbviol = sum(ineqv<ineqLB)
+ ubviol = sum(ineqv>ineqUB)
+ if( lbviol > 0 | ubviol > 0 ){
+ vidx = c(which(ineqv<ineqLB), which(ineqv>ineqUB))
+ vidx = unique(vidx)
+ gosolnp_rndpars = gosolnp_rndpars[-c(vidx),,drop=FALSE]
+ lvx = length(vidx)
+ } else{
+ vidx = 0
+ lvx = 0
+ }
+ } else{
+ ineqv = t(apply(gosolnp_rndpars, 1, FUN = function(x){
+ names(x) = gosolnp_parnames
+ ineqfun(x, ...)} ))
+
+ # check lower and upper violations
+ lbviol = apply(ineqv, 1, FUN = function(x) sum(any(x<ineqLB)))
+ ubviol = apply(ineqv, 1, FUN = function(x) sum(any(x>ineqUB)))
+ if( any(lbviol > 0) | any(ubviol > 0) ){
+ vidx = c(which(lbviol>0), which(ubviol>0))
+ vidx = unique(vidx)
+ gosolnp_rndpars = gosolnp_rndpars[-c(vidx),,drop=FALSE]
+ lvx = length(vidx)
+
+ } else{
+ vidx = 0
+ lvx = 0
+ }
+ }
+ if( trace ) cat(paste("\n...Excluded ", lvx, "/",n.restarts*n.sim, " Random Sequences\n", sep = ""))
+ }
+ # evaluate function value
+ if( trace ) cat("\nEvaluating Objective Function with Random Sampled Parameters...")
+ if( !is.null(cluster) ){
+ nx = dim(gosolnp_rndpars)[1]
+ clusterExport(cluster, c("gosolnp_rndpars", ".safefun"), envir = environment())
+ evfun = parLapply(cluster, as.list(1:nx), fun = function(i){
+ .safefun(gosolnp_rndpars[i, ], fun, gosolnp_parnames, ...)
+ })
+ evfun = as.numeric( unlist(evfun) )
+ } else{
+ evfun = apply(gosolnp_rndpars, 1, FUN = function(x) .safefun(x, fun, gosolnp_parnames, ...))
+ }
+ if( trace ) cat("ok!\n")
+ if( trace ) cat("\nSorting and Choosing Best Candidates for starting Solver...")
+ z = sort.int(evfun, index.return = T)
+ ans = gosolnp_rndpars[z$ix[1:n.restarts],,drop = FALSE]
+ prtable = cbind(ans, z$x[1:n.restarts])
+ if( trace ) cat("ok!\n")
+ colnames(prtable) = c(paste("par", 1:N, sep = ""), "objf")
+ if( trace ){
+ cat("\nStarting Parameters and Starting Objective Function:\n")
+ if(n.restarts == 1) print(t(prtable), digits = 4) else print(prtable, digits = 4)
+ }
+ return(prtable)
+}
+
+# form a barrier function before passing the parameters
+.randpars2 = function(pars, fixed, fun, eqfun, eqB, ineqfun, ineqLB, ineqUB, LB,
+ UB, distr, distr.opt, n.restarts, n.sim, rseed, trace = TRUE,
+ gosolnp_parnames, cluster, ...)
+{
+ if( trace ) cat("\nCalculating Random Initialization Parameters...")
+ N = length(LB)
+ gosolnp_idx = "a"
+ gosolnp_R = NULL
+ if(!is.null(ineqfun) && is.null(eqfun) ){
+ gosolnp_idx = "b"
+ gosolnp_R = 100
+ }
+ if( is.null(ineqfun) && !is.null(eqfun) ){
+ gosolnp_idx = "c"
+ gosolnp_R = 100
+ }
+ if(!is.null(ineqfun) && !is.null(eqfun) ){
+ gosolnp_idx = "d"
+ gosolnp_R = c(100,100)
+ }
+ gosolnp_rndpars = matrix(NA, ncol = N, nrow = n.sim * n.restarts)
+ if(!is.null(fixed)) for(i in 1:length(fixed)) gosolnp_rndpars[,fixed[i]] = pars[fixed[i]]
+ nf = 1:N
+ if(!is.null(fixed)) nf = nf[-c(fixed)]
+ gosolnp_m = length(nf)
+ set.seed(rseed)
+ for(i in 1:gosolnp_m){
+ j = nf[i]
+ gosolnp_rndpars[,j] = switch(distr[j],
+ .distr1(LB[j], UB[j], n.restarts*n.sim),
+ .distr2(LB[j], UB[j], n.restarts*n.sim, mean = distr.opt[[j]]$mean, sd = distr.opt[[j]]$sd),
+ .distr3(n.restarts*n.sim, mean = distr.opt[[j]]$mean, sd = distr.opt[[j]]$sd)
+ )
+ }
+ if( trace ) cat("ok!\n")
+ # Barrier Function
+ pclfn = function(x){
+ z=x
+ z[x<=0] = 0
+ z[x>0] = (0.9+z[x>0])^2
+ z
+ }
+ .lagrfun = function(pars, m, idx, fun, eqfun = NULL, eqB = 0, ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, ...)
+ {
+ fn = switch(idx,
+ "a" = fun(pars[1:m], ...),
+ "b" = fun(pars[1:m], ...) + pars[m+1]* sum( pclfn( c(ineqLB - ineqfun(pars[1:m], ...), ineqfun(pars[1:m], ...) - ineqUB) ) ),
+ "c" = fun(pars[1:m], ...) + sum( (eqfun(pars[1:m], ...) - eqB )^2 / pars[m+1]),
+ "d" = fun(pars[1:m], ...) + sum( (eqfun(pars[1:m], ...) - eqB )^2 / pars[m+1]) + pars[m+2]* sum( pclfn( c(ineqLB - ineqfun(pars[1:m], ...), ineqfun(pars[1:m], ...) - ineqUB) ) ) )
+ return(fn)
+ }
+
+ # evaluate function value
+ if( trace ) cat("\nEvaluating Objective Function with Random Sampled Parameters...")
+ if( !is.null(cluster) ){
+ nx = dim(gosolnp_rndpars)[1]
+ clusterExport(cluster, c("gosolnp_rndpars", "gosolnp_m", "gosolnp_idx",
+ "gosolnp_R"), envir = environment())
+ clusterExport(cluster, c("pclfn", ".lagrfun"), envir = environment())
+ evfun = parallel::parLapply(cluster, as.list(1:nx), fun = function(i){
+ .lagrfun(c(gosolnp_rndpars[i,], gosolnp_R), gosolnp_m,
+ gosolnp_idx, fun, eqfun, eqB, ineqfun, ineqLB,
+ ineqUB, ...)
+ })
+ evfun = as.numeric( unlist(evfun) )
+ } else{
+ evfun = apply(gosolnp_rndpars, 1, FUN = function(x){
+ .lagrfun(c(x,gosolnp_R), gosolnp_m, gosolnp_idx, fun, eqfun,
+ eqB, ineqfun, ineqLB, ineqUB, ...)})
+ }
+ if( trace ) cat("ok!\n")
+ if( trace ) cat("\nSorting and Choosing Best Candidates for starting Solver...")
+ z = sort.int(evfun, index.return = T)
+ #distmat = dist(evfun, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
+ ans = gosolnp_rndpars[z$ix[1:n.restarts],,drop = FALSE]
+ prtable = cbind(ans, z$x[1:n.restarts])
+ colnames(prtable) = c(paste("par", 1:N, sep = ""), "objf")
+ if( trace ){
+ cat("\nStarting Parameters and Starting Objective Function:\n")
+ if(n.restarts == 1) print(t(prtable), digits = 4) else print(prtable, digits = 4)
+ }
+ return(prtable)
+}
+
+
+.distr1 = function(LB, UB, n)
+{
+ runif(n, min = LB, max = UB)
+}
+
+.distr2 = function(LB, UB, n, mean, sd)
+{
+ rtruncnorm(n, a = as.double(LB), b = as.double(UB), mean = as.double(mean), sd = as.double(sd))
+}
+
+.distr3 = function(n, mean, sd)
+{
+ rnorm(n, mean = mean, sd = sd)
+}
+
+.safefun = function(pars, fun, gosolnp_parnames, ...){
+ # gosolnp_parnames = get("gosolnp_parnames", envir = .env)
+ names(pars) = gosolnp_parnames
+ v = fun(pars, ...)
+ if(is.na(v) | !is.finite(v) | is.nan(v)) {
+ warning(paste("\ngosolnp-->warning: ", v , " detected in function call...check your function\n", sep = ""), immediate. = FALSE)
+ v = 1e24
+ }
+ v
+}
diff --git a/R/helpers.R b/R/helpers.R
new file mode 100644
index 0000000..9b245b9
--- /dev/null
+++ b/R/helpers.R
@@ -0,0 +1,351 @@
+#################################################################################
+##
+## R package Rsolnp by Alexios Ghalanos and Stefan Theussl Copyright (C) 2009-2013
+## This file is part of the R package Rsolnp.
+##
+## The R package Rsolnp is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## The R package Rsolnp is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+#################################################################################
+
+.eps = .Machine$double.eps
+
+.subnpmsg = function(m){
+ g1 = c("solnp-->")
+ m1 = paste("\n",g1, "Redundant constraints were found. Poor\n",
+ g1, "intermediate results may result.Suggest that you\n",
+ g1, "remove redundant constraints and re-OPTIMIZE\n", sep = "")
+ m2 = paste("\n",g1, "The linearized problem has no feasible\n",
+ g1, "solution. The problem may not be feasible.\n", sep = "")
+ m3 = paste("\n",g1, "Minor optimization routine did not converge in the \n",
+ g1, "specified number of minor iterations. You may need\n",
+ g1, "to increase the number of minor iterations. \n", sep = "")
+ ans = switch(m,
+ m1 = m1,
+ m2 = m2,
+ m3 = m3)
+ cat(ans)
+}
+
+
+
+.checkpars = function(pars, LB, UB, .env)
+{
+ if(is.null(pars))
+ stop("\nsolnp-->error: must supply starting parameters\n", call. = FALSE)
+ if(!is.null(LB)){
+ if(length(pars) != length(LB))
+ stop("\nsolnp-->error: LB length not equal to parameter length\n", call. = FALSE)
+ if(is.null(UB)) UB = rep(.Machine$double.xmax/2, length(LB))
+ } else{
+ LB = NULL
+ }
+ if(!is.null(UB)){
+ if(length(pars) != length(UB))
+ stop("\nsolnp-->error: UB length not equal to parameter length\n", call. = FALSE)
+ if(is.null(LB)) LB = rep(-.Machine$double.xmax/2, length(UB))
+ } else{
+ UB = NULL
+ }
+ if(!is.null(UB) && any(LB > UB))
+ stop("\nsolnp-->error: UB must be greater than LB\n", call. = FALSE)
+
+ if(!is.null(UB) && any(LB == UB))
+ warning("\nsolnp-->warning: Equal Lower/Upper Bounds Found. Consider\n
+ excluding fixed parameters.\n", call. = FALSE)
+ # deal with infinite values as these are not accepted by solve.QP
+ if(!is.null(LB) && !any(is.finite(LB))){
+ idx = which(!is.finite(LB))
+ LB[idx] = sign(LB[idx])*.Machine$double.xmax/2
+ }
+ if(!is.null(UB) && !any(is.finite(UB))){
+ idx = which(!is.finite(UB))
+ UB[idx] = sign(UB[idx])*.Machine$double.xmax/2
+ }
+ assign(".LB", LB, envir = .env)
+ assign(".UB", UB, envir = .env)
+ return(1)
+}
+
+.checkfun = function(pars, fun, .env, ...)
+{
+ if(!is.function(fun)) stop("\nsolnp-->error: fun does not appear to be a function\n", call. = FALSE)
+ val = fun(pars, ...)
+ if(length(val) != 1) stop("\nsolnp-->error: objective function returns value of length greater than 1!\n", call. = FALSE)
+
+ assign(".solnp_fun", fun, envir = .env)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+ return(val)
+}
+
+# Might eventually use this, but really the user must take care of such problems
+# in their own function/setup
+.safefunx = function(pars, fun, .env, ...){
+ xnames = get("xnames", envir = .env)
+ names(pars) = xnames
+ v = fun(pars, ...)
+ if(is.na(v) | !is.finite(v) | is.nan(v)) {
+ warning(paste("\nsolnp-->warning: ", v , " detected in function call...check your function\n", sep = ""), immediate. = FALSE)
+ v = 1e24
+ }
+ v
+}
+
+.checkgrad = function(pars, fun, .env, ...)
+{
+ n = length(pars)
+ val = fun(pars, ...)
+ if(length(val)!=n)
+ stop("\nsolnp-->error: gradient vector length must be equal to length(pars)\n", call. = FALSE)
+ assign(".solnp_gradfun", fun, envir = .env)
+ return(val)
+}
+
+.checkhess = function(pars, fun, .env, ...)
+{
+ n = length(pars)
+ val = fun(pars, ...)
+ if(length(as.vector(val)) != (n*n))
+ stop("\nsolnp-->error: hessian must be of length length(pars) x length(pars)\n", call. = FALSE)
+ assign(".solnp_hessfun", fun, envir = .env)
+ return(val)
+}
+
+.checkineq = function(pars, fun, ineqLB, ineqUB, .env, ...)
+{
+ xnames = get("xnames", envir = .env)
+ val = fun(pars, ...)
+ n = length(val)
+ if(!is.null(ineqLB)){
+ if(length(ineqLB) != n)
+ stop("\nsolnp-->error: inequality function returns vector of different length to
+ inequality lower bounds\n", call. = FALSE)
+ } else{
+ stop("\nsolnp-->error: inequality function given without lower bounds\n", call. = FALSE)
+ }
+ if(!is.null(ineqUB)){
+ if(length(ineqUB) != n)
+ stop("\nsolnp-->error: inequality function returns vector of different length to
+ inequality upper bounds\n", call. = FALSE)
+ } else{
+ stop("\nsolnp-->error: inequality function given without upper bounds\n", call. = FALSE)
+ }
+ if(any(ineqLB > ineqUB))
+ stop("\nsolnp-->error: ineqUB must be greater than ineqLB\n", call. = FALSE)
+
+ assign(".ineqLB", ineqLB, envir = .env)
+ assign(".ineqUB", ineqUB, envir = .env)
+ .solnp_ineqfun = function(x, ...){
+ names(x) = xnames
+ fun(x, ...)
+ }
+ assign(".solnp_ineqfun", .solnp_ineqfun, envir = .env)
+ return(val)
+}
+
+
+.checkeq = function(pars, fun, eqB, .env, ...)
+{
+ xnames = get("xnames", envir = .env)
+ n = length(eqB)
+ val = fun(pars, ...) - eqB
+ if(length(val)!=n)
+ stop("\nsolnp-->error: equality function returns vector of different length
+ to equality value\n", call. = FALSE)
+ .eqB = eqB
+ assign(".eqB", .eqB, envir = .env)
+ .solnp_eqfun = function(x, ...){
+ names(x) = xnames
+ fun(x, ...) - .eqB
+ }
+ assign(".solnp_eqfun", .solnp_eqfun, envir = .env)
+ return(val)
+}
+
+
+# check the jacobian of inequality
+.cheqjacineq = function(pars, fun, .env, ...)
+{
+ # must be a matrix -> nrows = no.inequalities, ncol = length(pars)
+ val = fun(pars, ...)
+ .ineqLB = get(".ineqLB", envir = .env)
+ .ineqUB = get(".ineqUB", envir = .env)
+ if(!is.matrix(val))
+ stop("\nsolnp-->error: Jacobian of Inequality must return a matrix type object\n", call. = FALSE)
+ nd = dim(val)
+ if(nd[2] != length(pars))
+ stop("\nsolnp-->error: Jacobian of Inequality column dimension must be equal to length
+ of parameters\n", call. = FALSE)
+ if(nd[1] != length(.ineqUB))
+ stop("\nsolnp-->error: Jacobian of Inequality row dimension must be equal to length
+ of inequality bounds vector\n", call. = FALSE)
+ # as in inequality function, transforms from a 2 sided inequality to a one sided inequality
+ # (for the jacobian).
+ .solnp_ineqjac = function(x, ...) { retval = fun(x, ...); rbind( retval ) }
+ assign(".solnp_ineqjac", .solnp_ineqjac, envir = .env)
+ return(val)
+}
+
+# check the jacobian of equality
+.cheqjaceq = function(pars, fun, .env, ...)
+{
+ # must be a matrix -> nrows = no.equalities, ncol = length(pars)
+ val = fun(pars, ...)
+ .eqB = get(".eqB", envir = .env)
+ if(!is.matrix(val))
+ stop("\nsolnp-->error: Jacobian of Equality must return a matrix type object\n", call. = FALSE)
+ nd = dim(val)
+ if(nd[2] != length(pars))
+ stop("\nsolnp-->error: Jacobian of Equality column dimension must be equal to length of parameters\n", call. = FALSE)
+ if(nd[1] != length(.eqB))
+ stop("\nsolnp-->error: Jacobian of Equality row dimension must be equal to length of equality bounds vector\n", call. = FALSE)
+ assign(".solnp_eqjac", fun, envir = .env)
+ return(val)
+}
+
+# reporting function
+.report = function(iter, funv, pars)
+{
+ cat( paste( "\nIter: ", iter ," fn: ", format(funv, digits = 4, scientific = 5, nsmall = 4, zero.print = TRUE), "\t Pars: ", sep=""),
+ format(pars, digits = 4, scientific = 6, nsmall = 5, zero.print = TRUE) )
+}
+
+# finite difference gradient
+.fdgrad = function(pars, fun, ...)
+{
+ if(!is.null(fun)){
+
+ y0 = fun(pars, ...)
+ nx = length(pars)
+ grd = rep(0, nx)
+ deltax = sqrt(.eps)
+ for(i in 1:nx)
+ {
+ init = pars[i]
+ pars[i]= pars[i] + deltax
+ grd[i] = (fun(pars, ...) - y0) / deltax
+ pars[i] = init
+ }
+ }
+ else
+ {
+ grd = 0
+ }
+ return(grd)
+}
+
+# finite difference jacobian
+.fdjac = function(pars, fun, ...)
+{
+ nx = length(pars)
+ if(!is.null(fun))
+ {
+ y0 = fun(pars, ...)
+ nf = length (y0)
+ jac = matrix(0, nrow = nf, ncol= nx)
+ deltax = sqrt (.eps)
+ for(i in 1:nx)
+ {
+ init = pars[i]
+ pars[i]= pars[i] + deltax
+ jac[,i] = (fun(pars, ...) - y0) / deltax
+ pars[i] = init
+ }
+ } else{
+ jac = rep(0, nx)
+ }
+ return(jac)
+}
+
+.emptygrad = function(pars, ...)
+{
+ matrix(0, nrow = 0, ncol = 1)
+}
+
+.emptyjac = function(pars, ...)
+{
+ #matrix(0, nrow = 0, ncol = length(pars))
+ NULL
+}
+
+.emptyfun = function(pars, ...)
+{
+ NULL
+}
+
+.ineqlbfun = function(pars, .env, ...)
+{
+ LB = get(".solnp_LB", envir = .env)
+ UB = get(".solnp_UB", envir = .env)
+ .solnp_ineqfun = get(".solnp_ineqfun", envir = .env)
+ res = c(pars - LB, UB - pars)
+ if(!is.null(.solnp_ineqfun)) res = c(.solnp_ineqfun(pars, ...), res)
+ res
+}
+
+.ineqlbjac = function(pars, .env, ...)
+{
+ .solnp_ineqjac = get(".solnp_ineqjac", envir = .env)
+ n = length(pars)
+ res = rbind(diag(n), -diag(n))
+ if(!is.null(.solnp_ineqjac)) res = rbind(.solnp_ineqjac(pars, ...), res)
+ res
+}
+
+.solnpctrl = function(control){
+ # parameters check is now case independent
+ ans = list()
+ params = unlist(control)
+ if(is.null(params)) {
+ ans$rho = 1
+ ans$outer.iter = 400
+ ans$inner.iter = 800
+ ans$delta = 1.0e-7
+ ans$tol = 1.0e-8
+ ans$trace = 1
+ } else{
+ npar = tolower(names(unlist(control)))
+ names(params) = npar
+ if(any(substr(npar, 1, 3) == "rho")) ans$rho = as.numeric(params["rho"]) else ans$rho = 1
+ if(any(substr(npar, 1, 10) == "outer.iter")) ans$outer.iter = as.numeric(params["outer.iter"]) else ans$outer.iter = 400
+ if(any(substr(npar, 1, 10) == "inner.iter")) ans$inner.iter = as.numeric(params["inner.iter"]) else ans$inner.iter = 800
+ if(any(substr(npar, 1, 5) == "delta")) ans$delta = as.numeric(params["delta"]) else ans$delta = 1.0e-7
+ if(any(substr(npar, 1, 3) == "tol")) ans$tol = as.numeric(params["tol"]) else ans$tol = 1.0e-8
+ if(any(substr(npar, 1, 5) == "trace")) ans$trace = as.numeric(params["trace"]) else ans$trace = 1
+ }
+ return(ans)
+}
+
+.zeros = function( n = 1, m = 1)
+{
+ if(missing(m)) m = n
+ sol = matrix(0, nrow = n, ncol = m)
+ return(sol)
+}
+
+.ones = function(n = 1, m = 1)
+{
+ if(missing(m)) m = n
+ sol = matrix(1, nrow = n, ncol = m)
+ return(sol)
+}
+
+.vnorm = function(x)
+{
+ sum((x)^2)^(1/2)
+}
+
+.solvecond = function(x)
+{
+ z = svd(x)$d
+ if(any( z == 0 )) ret = Inf else ret = max( z ) / min( z )
+ return(ret)
+}
diff --git a/R/solnp.R b/R/solnp.R
new file mode 100644
index 0000000..4c6dea8
--- /dev/null
+++ b/R/solnp.R
@@ -0,0 +1,349 @@
+#################################################################################
+##
+## R package Rsolnp by Alexios Ghalanos and Stefan Theussl Copyright (C) 2009-2013
+## This file is part of the R package Rsolnp.
+##
+## The R package Rsolnp is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## The R package Rsolnp is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+#################################################################################
+
+# Based on the original solnp by Yinyu Ye
+# http://www.stanford.edu/~yyye/Col.html
+
+
+#----------------------------------------------------------------------------------
+# The Function SOLNP solves nonlinear programs in standard form:
+#
+# minimize J(P)
+# subject to EC(P) =0
+# IB(:,1)<= IC(P) <=IB(:,2)
+# PB(:,1)<= P <=PB(:,2).
+#where
+#
+# J : Cost objective scalar function
+# EC : Equality constraint vector function
+# IC : Inequality constraint vector function
+# P : Decision parameter vector
+# IB, PB : lower and upper bounds for IC and P.
+#----------------------------------------------------------------------------------
+
+# control list
+# RHO : penalty parameter
+# MAJIT: maximum number of major iterations
+# MINIT: maximum number of minor iterations
+# DELTA: relative step size in forward difference evaluation
+# TOL : tolerance on feasibility and optimality
+# defaults RHO=1, MAJIT=10, MINIT=10, DELTA=1.0e-5, TOL=1.0e-4
+
+solnp = function(pars, fun, eqfun = NULL, eqB = NULL, ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, LB = NULL, UB = NULL, control = list(), ...)
+{
+ # start timer
+ tic = Sys.time()
+ xnames = names(pars)
+ # get environment
+ .solnpenv <- environment()
+ assign("xnames", xnames, envir = .solnpenv)
+ # initiate function count
+ assign(".solnp_nfn", 0, envir = .solnpenv)
+ assign(".solnp_errors", 0, envir = .solnpenv)
+
+ # index of function indicators
+ # [1] length of pars
+ # [2] has function gradient?
+ # [3] has hessian?
+ # [4] has ineq?
+ # [5] ineq length
+ # [6] has jacobian (inequality)
+ # [7] has eq?
+ # [8] eq length
+ # [9] has jacobian (equality)
+ # [10] has upper / lower bounds
+ # [11] has either lower/upper bounds or ineq
+
+
+ ind = rep(0, 11)
+ np = ind[1] = length(pars)
+ # lower parameter bounds - indicator
+ # lpb[1]=1 means lower/upper bounds present
+ # lpb[2]=1 means lower/upper bounds OR inequality bounds present
+
+ # do parameter and LB/UB checks
+ check1 = .checkpars(pars, LB, UB, .solnpenv)
+
+ # .LB and .UB assigned
+
+ .LB = get(".LB", envir = .solnpenv)
+ .UB = get(".UB", envir = .solnpenv)
+
+
+ if( !is.null(.LB) || !is.null(.UB) ) ind[10] = 1
+
+ # do function checks and return starting value
+ funv = .checkfun(pars, fun, .solnpenv, ...)
+ #.solnp_fun assigned
+ .solnp_fun = get(".solnp_fun", envir = .solnpenv)
+
+ # Analytical Gradient Functionality not yet implemented in subnp function
+
+ # gradient and hessian checks
+ #if(!is.null(grad)){
+ # gradv = .checkgrad(pars, grad, .solnpenv, ...)
+ # ind[2] = 1
+ #} else{
+ # .solnp_gradfun = function(pars, ...) .fdgrad(pars, fun = .solnp_fun, ...)
+ ind[2] = 0
+ # gradv = .solnp_gradfun(pars, ...)
+ #}
+ # .solnp_gradfun(pars, ...) assigned
+
+ .solnp_hessfun = NULL
+ ind[3] = 0
+ #hessv = NULL
+ # .solnp_hessfun(pars, ...) assigned
+
+ # do inequality checks and return starting values
+
+ if(!is.null(ineqfun)){
+ ineqv = .checkineq(pars, ineqfun, ineqLB, ineqUB, .solnpenv, ...)
+ ind[4] = 1
+ nineq = length(ineqLB)
+ ind[5] = nineq
+
+ # check for infinities/nans
+ .ineqLBx = .ineqLB
+ .ineqUBx = .ineqUB
+ .ineqLBx[!is.finite(.ineqLB)] = -1e10
+ .ineqUBx[!is.finite(.ineqUB)] = 1e10
+ ineqx0 = (.ineqLBx + .ineqUBx)/2
+ #if(!is.null(ineqgrad)){
+ # ineqjacv = .cheqjacineq(pars, gradineq, .ineqUB, .ineqLB, .solnpenv, ...)
+ # ind[6] = 1
+ #} else{
+ # .solnp_ineqjac = function(pars, ...) .fdjac(pars, fun = .solnp_ineqfun, ...)
+ ind[6] = 0
+ #ineqjacv = .solnp_ineqjac(pars, ...)
+ #}
+ } else{
+ .solnp_ineqfun = function(pars, ...) .emptyfun(pars, ...)
+ # .solnp_ineqjac = function(pars, ...) .emptyjac(pars, ...)
+ ineqv = NULL
+ ind[4] = 0
+ nineq = 0
+ ind[5] = 0
+ ind[6] = 0
+ ineqx0 = NULL
+ .ineqLB = NULL
+ .ineqUB = NULL
+ }
+ # .solnp_ineqfun and .solnp_ineqjac assigned
+ # .ineqLB and .ineqUB assigned
+ .solnp_ineqfun = get(".solnp_ineqfun", envir = .solnpenv)
+ .ineqLB = get(".ineqLB", envir = .solnpenv)
+ .ineqUB = get(".ineqUB", envir = .solnpenv)
+
+
+ # equality checks
+ if(!is.null(eqfun)){
+ eqv = .checkeq(pars, eqfun, eqB, .solnpenv, ...)
+ ind[7] = 1
+ .eqB = get(".eqB", envir = .solnpenv)
+ neq = length(.eqB)
+ ind[8] = neq
+ #if(!is.null(eqgrad)){
+ # eqjacv = .cheqjaceq(pars, gradeq, .solnpenv, ...)
+ # ind[9] = 1
+ #} else{
+ # .solnp_eqjac = function(pars, ...) .fdjac(pars, fun = .solnp_eqfun, ...)
+ # eqjacv = .solnp_eqjac(pars, ...)
+ ind[9] = 0
+ #}
+ } else {
+ eqv = NULL
+ #eqjacv = NULL
+ .solnp_eqfun = function(pars, ...) .emptyfun(pars, ...)
+ #.solnp_eqjac = function(pars, ...) .emptyjac(pars, ...)
+ ind[7] = 0
+ neq = 0
+ ind[8] = 0
+ ind[9] = 0
+ }
+ # .solnp_eqfun(pars, ...) and .solnp_eqjac(pars, ...) assigned
+ # .solnp_eqB assigned
+ .solnp_eqfun = get(".solnp_eqfun", envir = .solnpenv)
+
+ if( ind[ 10 ] || ind [ 4 ]) ind[ 11 ] = 1
+
+ # parameter bounds (pb)
+ pb = rbind( cbind(.ineqLB, .ineqUB), cbind(.LB, .UB) )
+
+ # check control list
+ ctrl = .solnpctrl( control )
+ rho = ctrl[[ 1 ]]
+ # maxit = outer iterations
+ maxit = ctrl[[ 2 ]]
+ # minit = inner iterations
+ minit = ctrl[[ 3 ]]
+ delta = ctrl[[ 4 ]]
+ tol = ctrl[[ 5 ]]
+ trace = ctrl[[ 6 ]]
+
+ # total constraints (tc) = no.inequality constraints + no.equality constraints
+ tc = nineq + neq
+
+ # initialize fn value and inequalities and set to NULL those not needed
+ j = jh = funv
+ tt = 0 * .ones(3, 1)
+
+ if( tc > 0 ) {
+ # lagrange multipliers (lambda)
+ lambda = 0 * .ones(tc, 1)
+ # constraint vector = [1:neq 1:nineq]
+ constraint = c(eqv, ineqv)
+ if( ind[4] ) {
+ tmpv = cbind(constraint[ (neq + 1):tc ] - .ineqLB, .ineqUB - constraint[ (neq + 1):tc ] )
+ testmin = apply( tmpv, 1, FUN = function( x ) min(x[ 1 ], x[ 2 ]) )
+ if( all(testmin > 0) ) ineqx0 = constraint[ (neq + 1):tc ]
+ constraint[ (neq + 1):tc ] = constraint[ (neq + 1):tc ] - ineqx0
+ }
+ tt[ 2 ] = .vnorm(constraint)
+ if( max(tt[ 2 ] - 10 * tol, nineq, na.rm = TRUE) <= 0 ) rho = 0
+ } else{
+ lambda = 0
+ }
+ # starting augmented parameter vector
+ p = c(ineqx0, pars)
+ hessv = diag(np + nineq)
+ mu = np
+ .solnp_iter = 0
+ ob = c(funv, eqv, ineqv)
+
+ while( .solnp_iter < maxit ){
+ .solnp_iter = .solnp_iter + 1
+ .subnp_ctrl = c(rho, minit, delta, tol, trace)
+
+ # make the scale for the cost, the equality constraints, the inequality
+ # constraints, and the parameters
+ if( ind[7] ) {
+ # [1 neq]
+ vscale = c( ob[ 1 ], rep(1, neq) * max( abs(ob[ 2:(neq + 1) ]) ) )
+ } else {
+ vscale = 1
+ }
+
+ if( !ind[ 11 ] ) {
+ vscale = c(vscale, p)
+ } else {
+ # [ 1 neq np]
+ vscale = c(vscale, rep( 1, length.out = length(p) ) )
+ }
+
+ vscale = apply( matrix(vscale, ncol = 1), 1, FUN = function( x ) min( max( abs(x), tol ), 1/tol ) )
+
+ res = .subnp(pars = p, yy = lambda, ob = ob, hessv = hessv, lambda = mu, vscale = vscale,
+ ctrl = .subnp_ctrl, .env = .solnpenv, ...)
+ if(get(".solnp_errors", envir = .solnpenv) == 1){
+ maxit = .solnp_iter
+ }
+ p = res$p
+ lambda = res$y
+ hessv = res$hessv
+ mu = res$lambda
+ temp = p[ (nineq + 1):(nineq + np) ]
+ funv = .safefunx(temp, .solnp_fun, .env = .solnpenv, ...)
+ ctmp = get(".solnp_nfn", envir = .solnpenv)
+ assign(".solnp_nfn", ctmp + 1, envir = .solnpenv)
+
+ tempdf = cbind(temp, funv)
+
+ if( trace ){
+ .report(.solnp_iter, funv, temp)
+ }
+
+ eqv = .solnp_eqfun(temp, ...)
+ ineqv = .solnp_ineqfun(temp, ...)
+
+ ob = c(funv, eqv, ineqv)
+
+ tt[ 1 ] = (j - ob[ 1 ]) / max(abs(ob[ 1 ]), 1)
+ j = ob[ 1 ]
+
+ if( tc > 0 ){
+ constraint = ob[ 2:(tc + 1) ]
+
+ if( ind[ 4 ] ){
+ tempv = rbind( constraint[ (neq + 1):tc ] - pb[ 1:nineq, 1 ],
+ pb[ 1:nineq, 2 ] - constraint[ (neq + 1):tc ] )
+
+ if( min(tempv) > 0 ) {
+ p[ 1:nineq ] = constraint[ (neq + 1):tc ]
+ }
+
+ constraint[ (neq + 1):tc ] = constraint[ (neq + 1):tc ] - p[ 1:nineq ]
+ }
+
+ tt[ 3 ] = .vnorm(constraint)
+
+ if( tt[ 3 ] < 10 * tol ) {
+ rho = 0
+ mu = min(mu, tol)
+ }
+
+ if( tt[ 3 ] < 5 * tt[ 2 ]) {
+ rho = rho/5
+ }
+
+ if( tt[ 3 ] > 10 * tt[ 2 ]) {
+ rho = 5 * max( rho, sqrt(tol) )
+ }
+
+ if( max( c( tol + tt[ 1 ], tt[ 2 ] - tt[ 3 ] ) ) <= 0 ) {
+ lambda = 0
+ hessv = diag( diag ( hessv ) )
+ }
+
+ tt[ 2 ] = tt[ 3 ]
+ }
+
+ if( .vnorm( c(tt[ 1 ], tt[ 2 ]) ) <= tol ) {
+ maxit = .solnp_iter
+ }
+
+ jh = c(jh, j)
+ }
+
+ if( ind[ 4 ] ) {
+ ineqx0 = p[ 1:nineq ]
+ }
+
+ p = p[ (nineq + 1):(nineq + np) ]
+
+ if(get(".solnp_errors", envir = .solnpenv) == 1){
+ convergence = 2
+ if( trace ) cat( paste( "\nsolnp--> Solution not reliable....Problem Inverting Hessian.\n", sep="" ) )
+ } else{
+ if( .vnorm( c(tt[ 1 ], tt[ 2 ]) ) <= tol ) {
+ convergence = 0
+ if( trace ) cat( paste( "\nsolnp--> Completed in ", .solnp_iter, " iterations\n", sep="" ) )
+ } else{
+ convergence = 1
+ if( trace ) cat( paste( "\nsolnp--> Exiting after maximum number of iterations\n",
+ "Tolerance not achieved\n", sep="" ) )
+ }
+ }
+ # end timer
+ ctmp = get(".solnp_nfn", envir = .solnpenv)
+ toc = Sys.time() - tic
+ names(p) = xnames
+ ans = list(pars = p, convergence = convergence, values = as.numeric(jh), lagrange = lambda,
+ hessian = hessv, ineqx0 = ineqx0, nfuneval = ctmp, outer.iter = .solnp_iter,
+ elapsed = toc, vscale = vscale)
+ return( ans )
+}
diff --git a/R/subnp.R b/R/subnp.R
new file mode 100644
index 0000000..b1538f6
--- /dev/null
+++ b/R/subnp.R
@@ -0,0 +1,563 @@
+#################################################################################
+##
+## R package Rsolnp by Alexios Ghalanos and Stefan Theussl Copyright (C) 2009-2013
+## This file is part of the R package Rsolnp.
+##
+## The R package Rsolnp is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## The R package Rsolnp is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+#################################################################################
+
+# Based on the original subnp Yinyu Ye
+# http://www.stanford.edu/~yyye/Col.html
+.subnp = function(pars, yy, ob, hessv, lambda, vscale, ctrl, .env, ...)
+{
+ .solnp_fun = get(".solnp_fun", envir = .env)
+ .solnp_eqfun = get(".solnp_eqfun", envir = .env)
+ .solnp_ineqfun = get(".solnp_ineqfun", envir = .env)
+ ineqLB = get(".ineqLB", envir = .env)
+ ineqUB = get(".ineqUB", envir = .env)
+ LB = get(".LB", envir = .env)
+ UB = get(".UB", envir = .env)
+ #.solnp_gradfun = get(".solnp_gradfun", envir = .env)
+ #.solnp_eqjac = get(".solnp_eqjac", envir = .env)
+ #.solnp_ineqjac = get(".solnp_ineqjac", envir = .env)
+ ind = get("ind", envir = .env)
+
+ # pars [nineq + np]
+ rho = ctrl[ 1 ]
+ maxit = ctrl[ 2 ]
+ delta = ctrl[ 3 ]
+ tol = ctrl[ 4 ]
+ trace = ctrl[ 5 ]
+ # [1] length of pars
+ # [2] has function gradient?
+ # [3] has hessian?
+ # [4] has ineq?
+ # [5] ineq length
+ # [6] has jacobian (inequality)
+ # [7] has eq?
+ # [8] eq length
+ # [9] has jacobian (equality)
+ # [10] has upper / lower bounds
+ # [11] has either lower/upper bounds or ineq
+
+
+ neq = ind[ 8 ]
+ nineq = ind[ 5 ]
+ np = ind[ 1 ]
+ ch = 1
+ alp = c(0,0,0)
+ nc = neq + nineq
+ npic = np + nineq
+ p0 = pars
+
+ # pb [ 2 x (nineq + np) ]
+ pb = rbind( cbind(ineqLB, ineqUB), cbind(LB,UB) )
+ sob = numeric()
+ ptt = matrix()
+ sc = numeric()
+
+ # scale the cost, the equality constraints, the inequality constraints,
+ # the parameters (inequality parameters AND actual parameters),
+ # and the parameter bounds if there are any
+ # Also make sure the parameters are no larger than (1-tol) times their bounds
+ # ob [ 1 neq nineq]
+
+ ob = ob / vscale[ 1:(nc + 1) ]
+ # p0 [np]
+ p0 = p0 / vscale[ (neq + 2):(nc + np + 1) ]
+ if( ind[ 11 ] ) {
+
+ if( !ind[ 10 ] ) {
+ mm = nineq
+ } else {
+ mm = npic
+ }
+
+ pb = pb / cbind(vscale[ (neq + 2):(neq + mm + 1) ], vscale[ (neq + 2):(neq + mm + 1) ])
+ }
+
+ # scale the lagrange multipliers and the Hessian
+
+ if( nc > 0) {
+ # yy [total constraints = nineq + neq]
+ # scale here is [tc] and dot multiplied by yy
+ yy = vscale[ 2:(nc + 1) ] * yy / vscale[ 1 ]
+ }
+ # yy = [zeros 3x1]
+
+ # h is [ (np+nineq) x (np+nineq) ]
+ #columnvector %*% row vector (size h) then dotproduct h then dotdivide scale[1]
+ hessv = hessv * (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1)]) ) / vscale[ 1 ]
+ # h[ 8x8 eye]
+ j = ob[ 1 ]
+
+ if( ind[4] ) {
+
+ if( !ind[7] ) {
+ # [nineq x (nineq+np) ]
+ a = cbind( -diag(nineq), matrix(0, ncol = np, nrow = nineq) )
+ } else {
+ # [ (neq+nineq) x (nineq+np)]
+ a = rbind( cbind( 0 * .ones(neq, nineq), matrix(0, ncol = np, nrow = neq) ),
+ cbind( -diag(nineq), matrix(0, ncol = np, nrow = nineq) ) )
+ }
+
+ }
+ if( ind[7] && !ind[4] ) {
+ a = .zeros(neq, np)
+ }
+
+ if( !ind[7] && !ind[4] ) {
+ a = .zeros(1, np)
+ }
+
+ # gradient
+ g= 0 * .ones(npic, 1)
+ p = p0 [ 1:npic ]
+ if( nc > 0 ) {
+ # [ nc ]
+ constraint = ob[ 2:(nc + 1) ]
+ # constraint [5 0 11 3x1]
+ # gradient routine
+ for( i in 1:np ) {
+ # scale the parameters (non ineq)
+ p0[ nineq + i ] = p0[ nineq + i ] + delta
+ tmpv = p0[ (nineq + 1):npic ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env, ...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+
+ ob = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+ g[ nineq + i ] = (ob[ 1 ] - j) / delta
+ a[ , nineq + i ] = (ob[ 2:(nc + 1) ] - constraint) / delta
+ p0[ nineq + i ] = p0[ nineq + i ] - delta
+ }
+
+ if( ind[4] ) {
+ constraint[ (neq + 1):(neq + nineq) ] = constraint[ (neq + 1):(neq + nineq) ] - p0[ 1:nineq ]
+ }
+
+ # solver messages
+ if( .solvecond(a) > 1 / .eps ) {
+ if( trace ) .subnpmsg( "m1" )
+ }
+
+ # a(matrix) x columnvector - columnvector
+ # b [nc,1]
+ b = a %*% p0 - constraint
+ ch = -1
+ alp[ 1 ] = tol - max( abs(constraint) )
+
+ if( alp[ 1 ] <= 0 ) {
+ ch = 1
+
+ if( !ind[11] ) {
+ # a %*% t(a) gives [nc x nc]
+ # t(a) %*% above gives [(np+nc) x 1]
+ p0 = p0 - t(a) %*% solve(a %*% t(a), constraint)
+ alp[ 1 ] = 1
+ }
+
+ }
+
+ if( alp[ 1 ] <= 0 ) {
+ # this expands p0 to [nc+np+1]
+ p0[ npic + 1 ] = 1
+ a = cbind(a, -constraint)
+ # cx is rowvector
+ cx = cbind(.zeros(1,npic), 1)
+ dx = .ones(npic + 1, 1)
+ go = 1
+ minit = 0
+
+ while( go >= tol ) {
+ minit = minit + 1
+ # gap [(nc + np) x 2]
+ gap = cbind(p0[ 1:mm ] - pb[ , 1 ], pb[ , 2 ] - p0[ 1:mm ] )
+ # this sorts every row
+ gap = t( apply( gap, 1, FUN=function( x ) sort(x) ) )
+ dx[ 1:mm ] = gap[ , 1 ]
+ # expand dx by 1
+ dx[ npic + 1 ] = p0[ npic + 1 ]
+
+ if( !ind[10] ) {
+ dx[ (mm + 1):npic ] = max( c(dx[ 1:mm ], 100) ) * .ones(npic - mm, 1)
+ }
+ # t( a %*% diag( as.numeric(dx) ) ) gives [(np+nc + 1 (or more) x nc]
+ # dx * t(cx) dot product of columnvectors
+ # qr.solve returns [nc x 1]
+
+ # TODO: Catch errors here
+ y = try( qr.solve( t( a %*% diag( as.numeric(dx) , length(dx), length(dx) ) ), dx * t(cx) ), silent = TRUE)
+ if(inherits(y, "try-error")){
+ p = p0 * vscale[ (neq + 2):(nc + np + 1) ] # unscale the parameter vector
+ if( nc > 0 ) {
+ y = 0 # unscale the lagrange multipliers
+ }
+ hessv = vscale[ 1 ] * hessv / (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1) ]) )
+ ans = list(p = p, y = y, hessv = hessv, lambda = lambda)
+ assign(".solnp_errors", 1, envir = .env)
+ return(ans)
+ }
+ v = dx * ( dx *(t(cx) - t(a) %*% y) )
+
+ if( v[ npic + 1 ] > 0 ) {
+ z = p0[ npic + 1 ] / v[ npic + 1 ]
+
+ for( i in 1:mm ) {
+
+ if( v[ i ] < 0 ) {
+ z = min(z, -(pb[ i, 2 ] - p0[ i ]) / v[ i ])
+ } else if( v[ i ] > 0 ) {
+ z = min( z, (p0[ i ] - pb[ i , 1 ]) / v[ i ])
+ }
+ }
+
+ if( z >= p0[ npic + 1 ] / v[ npic + 1 ] ) {
+ p0 = p0 - z * v
+ } else {
+ p0 = p0 - 0.9 * z * v
+ }
+ go = p0[ npic + 1 ]
+
+ if( minit >= 10 ) {
+ go = 0
+ }
+
+ } else {
+ go = 0
+ minit = 10
+ }
+
+ }
+
+ if( minit >= 10 ) {
+ if( trace ) .subnpmsg( "m2" )
+ }
+
+ a = matrix(a[ , 1:npic ], ncol = npic)
+ b = a %*% p0[ 1:npic ]
+ }
+
+ }
+
+ p = p0 [ 1:npic ]
+ y = 0
+
+ if( ch > 0 ) {
+
+ tmpv = p[ (nineq + 1):npic ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env,...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+ ob = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+ }
+
+ j = ob[ 1 ]
+
+ if( ind[4] ) {
+ ob[ (neq + 2):(nc + 1) ] = ob[ (neq + 2):(nc + 1) ] - p[ 1:nineq ]
+
+ }
+
+ if( nc > 0 ) {
+ ob[ 2:(nc + 1) ] = ob[ 2:(nc + 1) ] - a %*% p + b
+ j = ob[ 1 ] - t(yy) %*% matrix(ob[ 2:(nc + 1) ], ncol=1) + rho * .vnorm(ob[ 2:(nc + 1) ]) ^ 2
+ }
+
+ minit = 0
+ while( minit < maxit ) {
+ minit = minit + 1
+
+ if( ch > 0 ) {
+
+ for( i in 1:np ) {
+
+ p[ nineq + i ] = p[ nineq + i ] + delta
+ tmpv = p[ (nineq + 1):npic ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env, ...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+ obm = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+
+ if( ind[4] ) {
+ obm[ (neq + 2):(nc + 1)] = obm[ (neq + 2):(nc + 1) ] - p[ 1:nineq ]
+ }
+
+ if( nc > 0 ) {
+
+ obm[ 2:(nc + 1) ] = obm[ 2:(nc + 1) ] - a %*% p + b
+ obm = obm[ 1 ] - t(yy) %*% obm[ 2:(nc + 1) ] + rho * .vnorm(obm[ 2:(nc + 1 ) ]) ^ 2
+ }
+
+ g[ nineq + i ] = (obm - j) / delta
+ p[ nineq + i ] = p[ nineq + i ] - delta
+ }
+
+ if( ind[4] ) {
+ g[ 1:nineq ] = 0
+ }
+
+ }
+
+ if( minit > 1 ) {
+ yg = g - yg
+ sx = p - sx
+ sc[ 1 ] = t(sx) %*% hessv %*% sx
+ sc[ 2 ] = t(sx) %*% yg
+
+ if( (sc[ 1 ] * sc[ 2 ]) > 0 ) {
+ sx = hessv %*% sx
+ hessv = hessv - ( sx %*% t(sx) ) / sc[ 1 ] + ( yg %*% t(yg) ) / sc[ 2 ]
+ }
+
+ }
+
+ dx = 0.01 * .ones(npic, 1)
+ if( ind[11] ) {
+
+ gap = cbind(p[ 1:mm ] - pb[ , 1 ], pb[ , 2 ] - p[ 1:mm ])
+ gap = t( apply( gap, 1, FUN = function( x ) sort(x) ) )
+ gap = gap[ , 1 ] + sqrt(.eps) * .ones(mm, 1)
+ dx[ 1:mm, 1 ] = .ones(mm, 1) / gap
+ if( !ind[10] ){
+ dx[ (mm + 1):npic, 1 ] = min (c( dx[ 1:mm, 1 ], 0.01) ) * .ones(npic - mm, 1)
+ }
+
+ }
+
+ go = -1
+ lambda = lambda / 10
+ while( go <= 0 ) {
+ cz = try(chol( hessv + lambda * diag( as.numeric(dx * dx), length(dx), length(dx) ) ), silent = TRUE)
+ if(inherits(cz, "try-error")){
+ p = p * vscale[ (neq + 2):(nc + np + 1) ] # unscale the parameter vector
+ if( nc > 0 ) {
+ y = 0 # unscale the lagrange multipliers
+ }
+ hessv = vscale[ 1 ] * hessv / (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1) ]) )
+ ans = list(p = p, y = y, hessv = hessv, lambda = lambda)
+ assign(".solnp_errors", 1, envir = .env)
+ return(ans)
+ }
+ cz = try(solve(cz), silent = TRUE)
+ if(inherits(cz, "try-error")){
+ p = p * vscale[ (neq + 2):(nc + np + 1) ] # unscale the parameter vector
+ if( nc > 0 ) {
+ y = 0 # unscale the lagrange multipliers
+ }
+ hessv = vscale[ 1 ] * hessv / (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1) ]) )
+ ans = list(p = p, y = y, hessv = hessv, lambda = lambda)
+ assign(".solnp_errors", 1, envir = .env)
+ return(ans)
+ }
+ yg = t(cz) %*% g
+
+ if( nc == 0 ) {
+ u = -cz %*% yg
+ } else{
+ y = try( qr.solve(t(cz) %*% t(a), yg), silent = TRUE )
+ if(inherits(y, "try-error")){
+ p = p * vscale[ (neq + 2):(nc + np + 1) ] # unscale the parameter vector
+ if( nc > 0 ) {
+ # y = vscale[ 1 ] * y / vscale[ 2:(nc + 1) ] # unscale the lagrange multipliers
+ y = 0
+ }
+ hessv = vscale[ 1 ] * hessv / (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1) ]) )
+ ans = list(p = p, y = y, hessv = hessv, lambda = lambda)
+ assign(".solnp_errors", 1, envir = .env)
+ return(ans)
+ }
+
+ u = -cz %*% (yg - ( t(cz) %*% t(a) ) %*% y)
+ }
+
+ p0 = u[ 1:npic ] + p
+ if( !ind[ 11 ] ) {
+ go = 1
+ } else {
+ go = min( c(p0[ 1:mm ] - pb[ , 1 ], pb[ , 2 ] - p0[ 1:mm ]) )
+ lambda = 3 * lambda
+ }
+
+ }
+
+ alp[ 1 ] = 0
+ ob1 = ob
+ ob2 = ob1
+ sob[ 1 ] = j
+ sob[ 2 ] = j
+ ptt = cbind(p, p)
+ alp[ 3 ] = 1.0
+ ptt = cbind(ptt, p0)
+ tmpv = ptt[ (nineq + 1):npic, 3 ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env, ...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+
+ ob3 = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+ sob[ 3 ] = ob3[ 1 ]
+
+ if( ind[4] ) {
+ ob3[ (neq + 2):(nc + 1) ] = ob3[ (neq + 2):(nc + 1) ] - ptt[ 1:nineq, 3 ]
+ }
+
+ if( nc > 0 ) {
+ ob3[ 2:(nc + 1) ] = ob3[ 2:(nc + 1) ] - a %*% ptt[ , 3 ] + b
+ sob[ 3 ] = ob3[ 1 ] - t(yy) %*% ob3[ 2:(nc + 1) ] + rho * .vnorm(ob3[ 2:(nc + 1) ]) ^ 2
+ }
+
+ go = 1
+ while( go > tol ) {
+ alp[ 2 ] = (alp[ 1 ] + alp[ 3 ]) / 2
+ ptt[ , 2 ] = (1 - alp[ 2 ]) * p + alp[ 2 ] * p0
+ tmpv = ptt[ (nineq + 1):npic, 2 ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env, ...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+
+ ob2 = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+
+ sob[ 2 ] = ob2[ 1 ]
+
+ if( ind[4] ) {
+ ob2[ (neq + 2):(nc + 1) ] = ob2[ (neq + 2):(nc + 1) ] - ptt[ 1:nineq , 2 ]
+ }
+
+ if( nc > 0 ) {
+ ob2[ 2:(nc + 1) ] = ob2[ 2:(nc + 1) ] - a %*% ptt[ , 2 ] + b
+ sob[ 2 ] = ob2[ 1 ] - t(yy) %*% ob2[ 2:(nc + 1) ] + rho * .vnorm(ob2[ 2:(nc + 1) ]) ^ 2
+ }
+
+ obm = max(sob)
+
+ if( obm < j ) {
+ obn = min(sob)
+ go = tol * (obm - obn) / (j - obm)
+ }
+
+ condif1 = sob[ 2 ] >= sob[ 1 ]
+ condif2 = sob[ 1 ] <= sob[ 3 ] && sob[ 2 ] < sob[ 1 ]
+ condif3 = sob[ 2 ] < sob[ 1 ] && sob[ 1 ] > sob[ 3 ]
+
+ if( condif1 ) {
+ sob[ 3 ] = sob[ 2 ]
+ ob3 = ob2
+ alp[ 3 ] = alp[ 2 ]
+ ptt[ , 3 ] = ptt[ , 2 ]
+ }
+
+ if( condif2 ) {
+ sob[ 3 ] = sob[ 2 ]
+ ob3 = ob2
+ alp[ 3 ] = alp[ 2 ]
+ ptt[ , 3 ] = ptt[ , 2 ]
+ }
+
+ if( condif3 ) {
+ sob[ 1 ] = sob[ 2 ]
+ ob1 = ob2
+ alp[ 1 ] = alp[ 2 ]
+ ptt[ , 1 ] = ptt[ , 2 ]
+ }
+
+ if( go >= tol ) {
+ go = alp[ 3 ] - alp[ 1 ]
+ }
+
+ }
+
+ sx = p
+ yg = g
+ ch = 1
+ obn = min(sob)
+ if( j <= obn ) {
+ maxit = minit
+ }
+
+ reduce = (j - obn) / ( 1 + abs(j) )
+
+ if( reduce < tol ) {
+ maxit = minit
+ }
+
+ condif1 = sob[ 1 ] < sob[ 2 ]
+ condif2 = sob[ 3 ] < sob[ 2 ] && sob[ 1 ] >= sob[ 2 ]
+ condif3 = sob[ 1 ] >= sob[ 2 ] && sob[ 3 ] >= sob[ 2 ]
+
+ if( condif1 ) {
+ j = sob[ 1 ]
+ p = ptt[ , 1 ]
+ ob = ob1
+ }
+
+ if( condif2 ) {
+ j = sob [ 3 ]
+ p = ptt[ , 3 ]
+ ob = ob3
+ }
+
+ if( condif3 ) {
+ j = sob[ 2 ]
+ p = ptt[ , 2 ]
+ ob = ob2
+ }
+
+ }
+
+ p = p * vscale[ (neq + 2):(nc + np + 1) ] # unscale the parameter vector
+
+ if( nc > 0 ) {
+ y = vscale[ 1 ] * y / vscale[ 2:(nc + 1) ] # unscale the lagrange multipliers
+ }
+
+ hessv = vscale[ 1 ] * hessv / (vscale[ (neq + 2):(nc + np + 1) ] %*% t(vscale[ (neq + 2):(nc + np + 1) ]) )
+ if( reduce > tol ) {
+ if( trace ) .subnpmsg( "m3" )
+ }
+
+ ans = list(p = p, y = y, hessv = hessv, lambda = lambda)
+ return( ans )
+}
+
+
+.fdgrad = function(i, p, delta, np, vscale, constraint, j, nineq, npic, nc, .solnp_fun,
+ .solnp_eqfun, .solnp_ineqfun, .env, ...)
+{
+ ans = list()
+ px = p
+ px[ nineq + i ] = px[ nineq + i ] + delta
+ tmpv = px[ (nineq + 1):npic ] * vscale[ (nc + 2):(nc + np + 1) ]
+ funv = .safefunx(tmpv, .solnp_fun, .env, ...)
+ eqv = .solnp_eqfun(tmpv, ...)
+ ineqv = .solnp_ineqfun(tmpv, ...)
+ ctmp = get(".solnp_nfn", envir = .env)
+ assign(".solnp_nfn", ctmp + 1, envir = .env)
+ ob = c(funv, eqv, ineqv) / vscale[ 1:(nc + 1) ]
+ ans$p = px
+ ans$ob = ob
+ ans$g = (ob[ 1 ] - j) / delta
+ ans$a = (ob[ 2:(nc + 1) ] - constraint) / delta
+ return( ans )
+}
+
+
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index cdb847e..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,14 +0,0 @@
-r-cran-rsolnp (1.16+dfsg-1) unstable; urgency=medium
-
- * New upstream version
- * xz compression in d/watch
- * Convert to dh-r
- * Canonical homepage for CRAN
-
- -- Andreas Tille <tille at debian.org> Sun, 06 Nov 2016 09:45:21 +0100
-
-r-cran-rsolnp (1.14+dfsg-1) unstable; urgency=low
-
- * Initial release (closes: #753125)
-
- -- Andreas Tille <tille at debian.org> Fri, 20 Jun 2014 11:47:43 +0200
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index ec63514..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-9
diff --git a/debian/control b/debian/control
deleted file mode 100644
index 052856e..0000000
--- a/debian/control
+++ /dev/null
@@ -1,24 +0,0 @@
-Source: r-cran-rsolnp
-Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
-Uploaders: Andreas Tille <tille at debian.org>
-Section: gnu-r
-Priority: optional
-Build-Depends: debhelper (>= 9),
- dh-r,
- r-base-dev,
- r-cran-truncnorm
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-cran-rsolnp/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-cran-rsolnp/trunk/
-Homepage: https://cran.r-project.org/package=Rsolnp
-
-Package: r-cran-rsolnp
-Architecture: any
-Depends: ${shlibs:Depends},
- ${misc:Depends},
- ${R:Depends}
-Recommends: ${R:Recommends}
-Suggests: ${R:Suggests}
-Description: GNU R general non-linear optimization
- This GNU R package provides general non-linear optimization using
- augmented lagrange multiplier method
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index 08c7d60..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,31 +0,0 @@
-Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: Rsolnp
-Upstream-Contact: Alexios Ghalanos <alexios at 4dscape.com>
-Source: https://cran.r-project.org/package=Rsolnp
-Files-Excluded: inst/doc/manual.ps
-
-Files: *
-Copyright: 2012-2016 Alexios Ghalanos, Stefan Theussl
-License: GPL-3+
-
-Files: debian/*
-Copyright: 2014-2016 Andreas Tille <tille at debian.org>
-License: GPL-3+
-
-License: GPL-3+
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- .
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
- .
- On Debian systems, the complete text of the GNU General Public
- License can be found in `/usr/share/common-licenses/GPL-3'.
-
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 68d9a36..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/usr/bin/make -f
-
-%:
- dh $@ --buildsystem R
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index a8f83d4..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,3 +0,0 @@
-version=3
-opts="repacksuffix=+dfsg,dversionmangle=s/\+dfsg//g,repack,compression=xz" \
- http://cran.r-project.org/src/contrib/Rsolnp_([-\d.]*)\.tar\.gz
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..3409145
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,32 @@
+citHeader("When using Rsolnp in publications, please cite both, the Rsolnp package, and the original algorithm:")
+
+## R >= 2.8.0 passes package metadata to citation().
+if(!exists("meta") || is.null(meta)) meta <- packageDescription("Rsolnp")
+
+year <- sub("-.*", "", meta$Date)
+note <- sprintf("R package version %s.", meta$Version)
+
+citEntry(entry="Manual",
+ title = "Rsolnp: General Non-linear Optimization Using Augmented
+ Lagrange Multiplier Method",
+ author = personList(as.person("Alexios Ghalanos and Stefan Theussl")),
+ year = year,
+ note = note,
+ textVersion =
+ paste("Alexios Ghalanos and Stefan Theussl",
+ sprintf("(%s).", year),
+ "Rsolnp: General Non-linear Optimization Using Augmented Lagrange Multiplier Method.",
+ note),
+ header = "To cite the Rsolnp package, please use:"
+)
+
+citEntry(entry="PhdThesis",
+ title = "Interior Algorithms for Linear, Quadratic, and Linearly Constrained Non-Linear Programming",
+ author = personList(as.person("Yinyu Ye")),
+ year = "1987",
+ school = "Department of {ESS}, Stanford University",
+ textVersion =
+ paste("Yinyu Ye (1987).",
+ "Interior Algorithms for Linear, Quadratic, and Linearly Constrained Non-Linear Programming. Ph.D. Thesis, Department of EES, Stanford University."),
+ header = "To cite the original SOLNP algorithm, please use:"
+)
diff --git a/inst/COPYRIGHTS b/inst/COPYRIGHTS
new file mode 100644
index 0000000..26f8fa3
--- /dev/null
+++ b/inst/COPYRIGHTS
@@ -0,0 +1,10 @@
+COPYRIGHT STATUS
+----------------
+
+This code is
+
+ Copyright (C) 2009-2013 by Alexios Ghalanos and Stefan Theussl.
+
+All code is subject to the GNU General Public License, Version 3. See
+the file COPYING for the exact conditions under which you may
+redistribute it.
diff --git a/inst/doc/manual.pdf b/inst/doc/manual.pdf
new file mode 100644
index 0000000..2fdfd5b
Binary files /dev/null and b/inst/doc/manual.pdf differ
diff --git a/man/Rsolnp-package.Rd b/man/Rsolnp-package.Rd
new file mode 100644
index 0000000..d3409c3
--- /dev/null
+++ b/man/Rsolnp-package.Rd
@@ -0,0 +1,25 @@
+\name{Rsolnp-package}
+\alias{Rsolnp-package}
+\alias{Rsolnp}
+\title{The Rsolnp package}
+\description{
+The Rsolnp package implements Y.Ye's general nonlinear augmented Lagrange
+multiplier method solver (SQP based solver).}
+\details{
+\tabular{ll}{
+Package: \tab Rsolnp\cr
+Type: \tab Package\cr
+Version: \tab 1.15\cr
+Date: \tab 2013-04-10\cr
+License: \tab GPL\cr
+LazyLoad: \tab yes\cr
+Depends: \tab stats, truncnorm, parallel\cr}
+}
+\author{
+Alexios Ghalanos and Stefan Theussl
+}
+\references{
+Y.Ye, \emph{Interior algorithms for linear, quadratic, and linearly constrained
+non linear programming}, PhD Thesis, Department of EES Stanford University,
+Stanford CA.
+}
\ No newline at end of file
diff --git a/man/benchmark.Rd b/man/benchmark.Rd
new file mode 100644
index 0000000..20f0716
--- /dev/null
+++ b/man/benchmark.Rd
@@ -0,0 +1,60 @@
+\name{benchmark}
+\Rdversion{1.1}
+\alias{benchmark}
+\title{
+The Rsolnp Benchmark Problems Suite.
+}
+\description{
+The function implements a set of benchmark problems against the MINOS solver of
+Murtagh and Saunders.
+}
+\usage{
+benchmark(id = "Powell")
+}
+\arguments{
+ \item{id}{
+The name of the benchmark problem. A call to the function \code{\link{benchmarkids}}
+will return the available benchmark problems.
+}
+}
+\details{
+The benchmarks were run on dual xeon server with 24GB of memory and windows 7
+operating system. The MINOS solver was used via the tomlab interface.
+}
+\value{
+A data.frame containing the benchmark data. The description of the benchmark
+problem can be accessed throught the \code{description} attribute of
+the data.frame.
+}
+\references{
+W.Hock and K.Schittkowski, \emph{Test Examples for Nonlinear Programming Codes},
+Lecture Notes in Economics and Mathematical Systems. Springer Verlag, 1981.\cr
+Y.Ye, \emph{Interior algorithms for linear, quadratic, and linearly constrained
+non linear programming}, PhD Thesis, Department of EES Stanford University,
+Stanford CA.\cr
+B.A.Murtagh and M.A.Saunders, \emph{MINOS 5.5 User's Guide, Report SOL 83-20R},
+Systems Optimization Laboratory, Stanford University (revised July 1998).\cr
+P. E. Gill, W. Murray, and M. A. Saunders, \emph{SNOPT An SQP algorithm for
+large-scale constrained optimization}, SIAM J. Optim., 12 (2002), pp.979-1006.
+}
+\author{
+Alexios Ghalanos and Stefan Theussl\cr
+Y.Ye (original matlab version of solnp)
+}
+\examples{
+\dontrun{
+benchmarkids()
+benchmark(id = "Powell")
+benchmark(id = "Alkylation")
+benchmark(id = "Box")
+benchmark(id = "RosenSuzuki")
+benchmark(id = "Wright4")
+benchmark(id = "Wright9")
+benchmark(id = "Electron")
+benchmark(id = "Permutation")
+# accessing the description
+test = benchmark(id = "Entropy")
+attr(test, "description")
+}
+}
+\keyword{optimize}
\ No newline at end of file
diff --git a/man/benchmarkids.Rd b/man/benchmarkids.Rd
new file mode 100644
index 0000000..5a32f78
--- /dev/null
+++ b/man/benchmarkids.Rd
@@ -0,0 +1,29 @@
+\name{benchmarkids}
+\alias{benchmarkids}
+\title{
+The Rsolnp Benchmark Problems Suite problem id's.
+}
+\description{
+Returns the id's of available benchmark in the Rsolnp Benchmark Problems Suite.
+}
+\usage{
+benchmarkids()
+}
+\value{
+A character vector of problem id's.
+}
+\references{
+W.Hock and K.Schittkowski, \emph{Test Examples for Nonlinear Programming Codes},
+Lecture Notes in Economics and Mathematical Systems. Springer Verlag, 1981.\cr
+Y.Ye, \emph{Interior algorithms for linear, quadratic, and linearly constrained
+non linear programming}, PhD Thesis, Department of EES Stanford University,
+Stanford CA.\cr
+}
+\author{
+Alexios Ghalanos and Stefan Theussl\cr
+Y.Ye (original matlab version of solnp)
+}
+\examples{
+benchmarkids()
+}
+\keyword{optimize}
\ No newline at end of file
diff --git a/man/gosolnp.Rd b/man/gosolnp.Rd
new file mode 100644
index 0000000..ebb41c2
--- /dev/null
+++ b/man/gosolnp.Rd
@@ -0,0 +1,247 @@
+\name{gosolnp}
+\alias{gosolnp}
+\title{
+Random Initialization and Multiple Restarts of the solnp solver.
+}
+\description{
+When the objective function is non-smooth or has many local minima, it is hard
+to judge the optimality of the solution, and this usually depends critically on
+the starting parameters. This function enables the generation of a set of
+randomly chosen parameters from which to initialize multiple restarts of the
+solver (see note for details).
+}
+\usage{
+gosolnp(pars = NULL, fixed = NULL, fun, eqfun = NULL, eqB = NULL, ineqfun = NULL,
+ineqLB = NULL, ineqUB = NULL, LB = NULL, UB = NULL, control = list(),
+distr = rep(1, length(LB)), distr.opt = list(), n.restarts = 1, n.sim = 20000,
+cluster = NULL, rseed = NULL, ...)
+}
+\arguments{
+ \item{pars}{
+The starting parameter vector. This is not required unless the fixed option is
+also used.
+}
+ \item{fixed}{
+The numeric index which indicates those parameters which should stay fixed
+instead of being randomly generated.
+}
+ \item{fun}{
+The main function which takes as first argument the parameter vector and returns
+a single value.
+}
+ \item{eqfun}{
+(Optional) The equality constraint function returning the vector of evaluated
+equality constraints.
+}
+ \item{eqB}{
+(Optional) The equality constraints.
+}
+ \item{ineqfun}{
+(Optional) The inequality constraint function returning the vector of evaluated
+inequality constraints.
+}
+ \item{ineqLB}{
+(Optional) The lower bound of the inequality constraints.
+}
+ \item{ineqUB}{
+(Optional) The upper bound of the inequality constraints.
+}
+ \item{LB}{
+The lower bound on the parameters. This is not optional in this function.
+}
+ \item{UB}{
+The upper bound on the parameters. This is not optional in this function.
+}
+ \item{control}{
+(Optional) The control list of optimization parameters. The \code{eval.type}
+option in this control list denotes whether to evaluate the function as is and
+exclude inequality violations in the final ranking (default, value = 1),
+else whether to evaluate a penalty barrier function comprised of the objective
+and all constraints (value = 2). See \code{solnp} function documentation for
+details of the remaining control options.
+}
+ \item{distr}{
+A numeric vector of length equal to the number of parameters, indicating the
+choice of distribution to use for the random parameter generation. Choices are
+uniform (1), truncated normal (2), and normal (3).
+}
+ \item{distr.opt}{
+If any choice in \code{distr} was anything other than uniform (1), this is a
+list equal to the length of the parameters with sub-components for the mean and
+sd, which are required in the truncated normal and normal distributions.
+}
+ \item{n.restarts}{
+The number of solver restarts required.
+}
+ \item{n.sim}{
+The number of random parameters to generate for every restart of the solver.
+Note that there will always be significant rejections if inequality bounds are
+present. Also, this choice should also be motivated by the width of the upper
+and lower bounds.
+}
+ \item{cluster}{
+If you want to make use of parallel functionality, initialize and pass a cluster
+object from the parallel package (see details), and remember to terminate it!
+}
+ \item{rseed}{
+(Optional) A seed to initiate the random number generator, else system time will
+be used.
+}
+ \item{\dots}{
+(Optional) Additional parameters passed to the main, equality or inequality
+functions
+}
+}
+\details{
+Given a set of lower and upper bounds, the function generates, for those
+parameters not set as fixed, random values from one of the 3 chosen
+distributions. Depending on the \code{eval.type} option of the \code{control}
+argument, the function is either directly evaluated for those points not
+violating any inequality constraints, or indirectly via a penalty barrier
+function jointly comprising the objective and constraints. The resulting values
+are then sorted, and the best N (N = random.restart) parameter vectors
+(corresponding to the best N objective function values) chosen in order to
+initialize the solver. Since version 1.14, it is up to the user to prepare and
+pass a cluster object from the parallel package for use with gosolnp, after
+which the parLapply function is used. If your function makes use of additional
+packages, or functions, then make sure to export them via the \code{clusterExport}
+function of the parallel package. Additional arguments passed to the solver via the
+\dots option are evaluated and exported by gosolnp to the cluster.
+}
+\value{
+A list containing the following values:
+\item{pars}{Optimal Parameters.}
+\item{convergence }{Indicates whether the solver has converged (0) or not (1).}
+\item{values}{Vector of function values during optimization with last one the
+value at the optimal.}
+\item{lagrange}{The vector of Lagrange multipliers.}
+\item{hessian}{The Hessian at the optimal solution.}
+\item{ineqx0}{The estimated optimal inequality vector of slack variables used
+for transforming the inequality into an equality constraint.}
+\item{nfuneval}{The number of function evaluations.}
+\item{elapsed}{Time taken to compute solution.}
+\item{start.pars}{The parameter vector used to start the solver}
+}
+\references{
+Y.Ye, \emph{Interior algorithms for linear, quadratic, and linearly constrained
+non linear programming}, PhD Thesis, Department of EES Stanford University,
+Stanford CA.\cr
+Hu, X. and Shonkwiler, R. and Spruill, M.C. \emph{Random Restarts in Global
+Optimization}, 1994, Georgia Institute of technology, Atlanta.
+}
+\author{
+Alexios Ghalanos and Stefan Theussl\cr
+Y.Ye (original matlab version of solnp)
+}
+\note{
+The choice of which distribution to use for randomly sampling the parameter
+space should be driven by the user's knowledge of the problem and confidence or
+lack thereof of the parameter distribution. The uniform distribution indicates
+a lack of confidence in the location or dispersion of the parameter, while the
+truncated normal indicates a more confident choice in both the location and
+dispersion. On the other hand, the normal indicates perhaps a lack of knowledge
+in the upper or lower bounds, but some confidence in the location and dispersion
+of the parameter. In using choices (2) and (3) for \code{distr},
+the \code{distr.opt} list must be supplied with \code{mean} and \code{sd} as
+subcomponents for those parameters not using the uniform (the examples section
+hopefully clarifies the usage).
+}
+\examples{
+\dontrun{
+# [Example 1]
+# Distributions of Electrons on a Sphere Problem:
+# Given n electrons, find the equilibrium state distribution (of minimal Coulomb
+# potential) of the electrons positioned on a conducting sphere. This model is
+# from the COPS benchmarking suite. See http://www-unix.mcs.anl.gov/~more/cops/.
+gofn = function(dat, n)
+{
+
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
+ jj = matrix(1:n, ncol = n, nrow = n)
+ ij = which(ii<jj, arr.ind = TRUE)
+ i = ij[,1]
+ j = ij[,2]
+ # Coulomb potential
+ potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
+ potential
+}
+
+goeqfn = function(dat, n)
+{
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ apply(cbind(x^2, y^2, z^2), 1, "sum")
+}
+
+n = 25
+LB = rep(-1, 3*n)
+UB = rep(1, 3*n)
+eqB = rep(1, n)
+ans = gosolnp(pars = NULL, fixed = NULL, fun = gofn, eqfun = goeqfn, eqB = eqB,
+LB = LB, UB = UB, control = list(outer.iter = 100, trace = 1),
+distr = rep(1, length(LB)), distr.opt = list(), n.restarts = 2, n.sim = 20000,
+rseed = 443, n = 25)
+# should get a function value around 243.813
+
+# [Example 2]
+# Parallel functionality for solving the Upper to Lower CVaR problem (not properly
+# formulated...for illustration purposes only).
+
+mu =c(1.607464e-04, 1.686867e-04, 3.057877e-04, 1.149289e-04, 7.956294e-05)
+sigma = c(0.02307198,0.02307127,0.01953382,0.02414608,0.02736053)
+R = matrix(c(1, 0.408, 0.356, 0.347, 0.378, 0.408, 1, 0.385, 0.565, 0.578, 0.356,
+0.385, 1, 0.315, 0.332, 0.347, 0.565, 0.315, 1, 0.662, 0.378, 0.578,
+0.332, 0.662, 1), 5,5, byrow=TRUE)
+# Generate Random deviates from the multivariate Student distribution
+set.seed(1101)
+v = sqrt(rchisq(10000, 5)/5)
+S = chol(R)
+S = matrix(rnorm(10000 * 5), 10000) \%*\% S
+ret = S/v
+RT = as.matrix(t(apply(ret, 1, FUN = function(x) x*sigma+mu)))
+# setup the functions
+.VaR = function(x, alpha = 0.05)
+{
+ VaR = quantile(x, probs = alpha, type = 1)
+ VaR
+}
+
+.CVaR = function(x, alpha = 0.05)
+{
+ VaR = .VaR(x, alpha)
+ X = as.vector(x[, 1])
+ CVaR = VaR - 0.5 * mean(((VaR-X) + abs(VaR-X))) / alpha
+ CVaR
+}
+.fn1 = function(x,ret)
+{
+ port=ret\%*\%x
+ obj=-.CVaR(-port)/.CVaR(port)
+ return(obj)
+}
+
+# abs(sum) of weights ==1
+.eqn1 = function(x,ret)
+{
+ sum(abs(x))
+}
+
+LB=rep(0,5)
+UB=rep(1,5)
+pars=rep(1/5,5)
+ctrl = list(delta = 1e-10, tol = 1e-8, trace = 0)
+cl = makePSOCKcluster(2)
+# export the auxilliary functions which are used and cannot be seen by gosolnp
+clusterExport(cl, c(".CVaR", ".VaR"))
+ans = gosolnp(pars, fun = .fn1, eqfun = .eqn1, eqB = 1, LB = LB, UB = UB,
+n.restarts = 2, n.sim=500, cluster = cl, ret = RT)
+ans
+# don't forget to stop the cluster!
+stopCluster(cl)
+}
+}
+\keyword{optimize}
diff --git a/man/solnp.Rd b/man/solnp.Rd
new file mode 100644
index 0000000..35eec33
--- /dev/null
+++ b/man/solnp.Rd
@@ -0,0 +1,131 @@
+\name{solnp}
+\alias{solnp}
+\title{
+Nonlinear optimization using augmented Lagrange method.
+}
+\description{
+The solnp function is based on the solver by Yinyu Ye which solves the general
+nonlinear programming problem:\cr
+\deqn{\min f(x)}{min f(x)}
+\deqn{\mathrm{s.t.}}{s.t.}
+\deqn{g(x) = 0}{g(x) = 0}
+\deqn{l_h \leq h(x) \leq u_h}{l[h] <= h(x) <= u[h]}
+\deqn{l_x \leq x \leq u_x}{l[x] <= x <= u[x]}
+where, \eqn{f(x)}, \eqn{g(x)} and \eqn{h(x)} are smooth functions.
+}
+\usage{
+solnp(pars, fun, eqfun = NULL, eqB = NULL, ineqfun = NULL, ineqLB = NULL,
+ineqUB = NULL, LB = NULL, UB = NULL, control = list(), ...)
+}
+\arguments{
+ \item{pars}{
+The starting parameter vector.
+}
+ \item{fun}{
+The main function which takes as first argument the parameter vector and returns
+a single value.
+}
+ \item{eqfun}{
+(Optional) The equality constraint function returning the vector of evaluated
+equality constraints.
+}
+ \item{eqB}{
+(Optional) The equality constraints.
+}
+ \item{ineqfun}{
+(Optional) The inequality constraint function returning the vector of evaluated
+inequality constraints.
+}
+ \item{ineqLB}{
+(Optional) The lower bound of the inequality constraints.
+}
+ \item{ineqUB}{
+(Optional) The upper bound of the inequality constraints.
+}
+ \item{LB}{
+(Optional) The lower bound on the parameters.
+}
+ \item{UB}{
+(Optional) The upper bound on the parameters.
+}
+ \item{control}{
+(Optional) The control list of optimization parameters. See below for details.
+}
+ \item{\dots}{
+(Optional) Additional parameters passed to the main, equality or inequality
+functions. Note that the main and constraint functions must take the exact same
+arguments, irrespective of whether they are used by all of them.
+}
+}
+\details{
+The solver belongs to the class of indirect solvers and implements the augmented
+Lagrange multiplier method with an SQP interior algorithm.
+}
+\value{
+A list containing the following values:
+\item{pars}{Optimal Parameters.}
+\item{convergence }{Indicates whether the solver has converged (0) or not
+(1 or 2).}
+\item{values}{Vector of function values during optimization with last one the
+value at the optimal.}
+\item{lagrange}{The vector of Lagrange multipliers.}
+\item{hessian}{The Hessian of the augmented problem at the optimal solution.}
+\item{ineqx0}{The estimated optimal inequality vector of slack variables used
+for transforming the inequality into an equality constraint.}
+\item{nfuneval}{The number of function evaluations.}
+\item{elapsed}{Time taken to compute solution.}
+}
+\section{Control}{
+\describe{
+\item{rho}{This is used as a penalty weighting scaler for infeasibility in the
+augmented objective function. The higher its value the more the weighting to
+bring the solution into the feasible region (default 1). However, very high
+values might lead to numerical ill conditioning or significantly slow down
+convergence.}
+\item{outer.iter}{Maximum number of major (outer) iterations (default 400).}
+\item{inner.iter}{Maximum number of minor (inner) iterations (default 800).}
+\item{delta}{Relative step size in forward difference evaluation
+(default 1.0e-7).}
+\item{tol}{ Relative tolerance on feasibility and optimality (default 1e-8).}
+\item{trace}{The value of the objective function and the parameters is printed
+at every major iteration (default 1).}
+}}
+\references{
+Y.Ye, \emph{Interior algorithms for linear, quadratic, and linearly constrained
+non linear programming}, PhD Thesis, Department of EES Stanford University,
+Stanford CA.
+}
+\author{
+Alexios Ghalanos and Stefan Theussl\cr
+Y.Ye (original matlab version of solnp)
+}
+\note{
+The control parameters \code{tol} and \code{delta} are key in getting any
+possibility of successful convergence, therefore it is suggested that the user
+change these appropriately to reflect their problem specification.\cr
+The solver is a local solver, therefore for problems with rough surfaces and
+many local minima there is absolutely no reason to expect anything other than a
+local solution.
+}
+\examples{
+# From the original paper by Y.Ye
+# see the unit tests for more....
+#---------------------------------------------------------------------------------
+# POWELL Problem
+fn1=function(x)
+{
+ exp(x[1]*x[2]*x[3]*x[4]*x[5])
+}
+
+eqn1=function(x){
+ z1=x[1]*x[1]+x[2]*x[2]+x[3]*x[3]+x[4]*x[4]+x[5]*x[5]
+ z2=x[2]*x[3]-5*x[4]*x[5]
+ z3=x[1]*x[1]*x[1]+x[2]*x[2]*x[2]
+ return(c(z1,z2,z3))
+}
+
+
+x0 = c(-2, 2, 2, -1, -1)
+powell=solnp(x0, fun = fn1, eqfun = eqn1, eqB = c(10, 0, -1))
+}
+\keyword{optimize}
\ No newline at end of file
diff --git a/man/startpars.Rd b/man/startpars.Rd
new file mode 100644
index 0000000..8dec0b8
--- /dev/null
+++ b/man/startpars.Rd
@@ -0,0 +1,172 @@
+\name{startpars}
+\alias{startpars}
+\title{
+Generates and returns a set of starting parameters by sampling the parameter
+space based on the evaluation of the function and constraints.
+}
+\description{
+A simple penalty barrier function is formed which is then evaluated at randomly
+sampled points based on the upper and lower parameter bounds
+(when \code{eval.type} = 2), else the objective function directly for values not
+violating any inequality constraints (when \code{eval.type} = 1). The sampled
+points can be generated from the uniform, normal or truncated normal
+distributions.
+}
+\usage{
+startpars(pars = NULL, fixed = NULL, fun, eqfun = NULL, eqB = NULL,
+ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, LB = NULL, UB = NULL,
+distr = rep(1, length(LB)), distr.opt = list(), n.sim = 20000, cluster = NULL,
+rseed = NULL, bestN = 15, eval.type = 1, trace = FALSE, ...)
+}
+\arguments{
+ \item{pars}{
+The starting parameter vector. This is not required unless the fixed option is
+also used.
+}
+ \item{fixed}{
+The numeric index which indicates those parameters which should stay fixed
+instead of being randomly generated.
+}
+ \item{fun}{
+The main function which takes as first argument the parameter vector and returns
+a single value.
+}
+ \item{eqfun}{
+(Optional) The equality constraint function returning the vector of evaluated
+equality constraints.
+}
+ \item{eqB}{
+(Optional) The equality constraints.
+}
+ \item{ineqfun}{
+(Optional) The inequality constraint function returning the vector of evaluated
+inequality constraints.
+}
+ \item{ineqLB}{
+(Optional) The lower bound of the inequality constraints.
+}
+ \item{ineqUB}{
+(Optional) The upper bound of the inequality constraints.
+}
+ \item{LB}{
+The lower bound on the parameters. This is not optional in this function.
+}
+ \item{UB}{
+The upper bound on the parameters. This is not optional in this function.
+}
+ \item{distr}{
+A numeric vector of length equal to the number of parameters, indicating the
+choice of distribution to use for the random parameter generation. Choices are
+uniform (1), truncated normal (2), and normal (3).
+}
+ \item{distr.opt}{
+If any choice in \code{distr} was anything other than uniform (1), this is a
+list equal to the length of the parameters with sub-components for the mean and
+sd, which are required in the truncated normal and normal distributions.
+}
+ \item{bestN}{
+The best N (less than or equal to n.sim) set of parameters to return.
+}
+ \item{n.sim}{
+The number of random parameter sets to generate.
+}
+ \item{cluster}{
+If you want to make use of parallel functionality, initialize and pass a cluster
+object from the parallel package (see details), and remember to terminate it!
+}
+ \item{rseed}{
+(Optional) A seed to initiate the random number generator, else system time will
+be used.
+}
+ \item{eval.type}{
+Either 1 (default) for the direction evaluation of the function (excluding
+inequality constraint violations) or 2 for the penalty barrier method.
+}
+ \item{trace}{
+(logical) Whether to display the progress of the function evaluation.
+}
+ \item{\dots}{
+(Optional) Additional parameters passed to the main, equality or inequality
+functions
+}
+}
+\details{
+Given a set of lower and upper bounds, the function generates, for those
+parameters not set as fixed, random values from one of the 3 chosen
+distributions. For simple functions with only inequality constraints, the direct
+method (\code{eval.type} = 1) might work better. For more complex setups with
+both equality and inequality constraints the penalty barrier method
+(\code{eval.type} = 2)might be a better choice.
+}
+\value{
+A matrix of dimension bestN x (no.parameters + 1). The last column is the
+evaluated function value.
+}
+\author{
+Alexios Ghalanos and Stefan Theussl\cr
+}
+\note{
+The choice of which distribution to use for randomly sampling the parameter
+space should be driven by the user's knowledge of the problem and confidence or
+lack thereof of the parameter distribution. The uniform distribution indicates a
+lack of confidence in the location or dispersion of the parameter, while the
+truncated normal indicates a more confident choice in both the location and
+dispersion. On the other hand, the normal indicates perhaps a lack
+of knowledge in the upper or lower bounds, but some confidence in the location
+and dispersion of the parameter. In using choices (2) and (3) for \code{distr},
+the \code{distr.opt} list must be supplied with \code{mean} and \code{sd} as
+subcomponents for those parameters not using the uniform.
+}
+\examples{
+\dontrun{
+library(Rsolnp)
+library(parallel)
+# Windows
+cl = makePSOCKcluster(2)
+# Linux:
+# makeForkCluster(nnodes = getOption("mc.cores", 2L), ...)
+
+gofn = function(dat, n)
+{
+
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
+ jj = matrix(1:n, ncol = n, nrow = n)
+ ij = which(ii<jj, arr.ind = TRUE)
+ i = ij[,1]
+ j = ij[,2]
+ # Coulomb potential
+ potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
+ potential
+}
+
+goeqfn = function(dat, n)
+{
+ x = dat[1:n]
+ y = dat[(n+1):(2*n)]
+ z = dat[(2*n+1):(3*n)]
+ apply(cbind(x^2, y^2, z^2), 1, "sum")
+}
+n = 25
+LB = rep(-1, 3*n)
+UB = rep( 1, 3*n)
+eqB = rep( 1, n)
+
+sp = startpars(pars = NULL, fixed = NULL, fun = gofn , eqfun = goeqfn,
+eqB = eqB, ineqfun = NULL, ineqLB = NULL, ineqUB = NULL, LB = LB, UB = UB,
+distr = rep(1, length(LB)), distr.opt = list(), n.sim = 2000,
+cluster = cl, rseed = 100, bestN = 15, eval.type = 2, n = 25)
+#stop cluster
+stopCluster(cl)
+# the last column is the value of the evaluated function (here it is the barrier
+# function since eval.type = 2)
+print(round(apply(sp, 2, "mean"), 3))
+# remember to remove the last column
+ans = solnp(pars=sp[1,-76],fun = gofn , eqfun = goeqfn , eqB = eqB, ineqfun = NULL,
+ineqLB = NULL, ineqUB = NULL, LB = LB, UB = UB, n = 25)
+# should get a value of around 243.8162
+}}
+
+\keyword{optimize}
\ No newline at end of file
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-rsolnp.git
More information about the debian-med-commit
mailing list