diff --git a/.Rproj.user/864BD013/cpp-compilation-config b/.Rproj.user/864BD013/cpp-compilation-config new file mode 100644 index 0000000..7cd788c --- /dev/null +++ b/.Rproj.user/864BD013/cpp-compilation-config @@ -0,0 +1,28 @@ +{ + "args": [ + "-std=gnu++11", + "-IC:/PROGRA~1/R/R-41~1.2/include", + "-DNDEBUG", + "-IC:/Users/vrfra/Documents/GitHub/optimg/inst/include", + "-DBOOST_NO_LONG_LONG", + "-DBOOST_NO_AUTO_PTR", + "-DBOOST_BIND_GLOBAL_PLACEHOLDERS", + "-DRCPP_NO_RTTI", + "-DRCPP_USE_UNWIND_PROTECT", + "-IC:/Users/vrfra/Documents/R/win-library/4.1/Rcpp/include", + "-IC:/Users/vrfra/Documents/R/win-library/4.1/RcppArmadillo/include", + "-nostdinc", + "-IC:/rtools40/mingw64/include/c++/8.3.0", + "-IC:/rtools40/mingw64/include/c++/8.3.0/x86_64-w64-mingw32", + "-IC:/rtools40/mingw64/include/c++/8.3.0/backward", + "-IC:/rtools40/mingw64/lib/gcc/x86_64-w64-mingw32/8.3.0/include", + "-IC:/rtools40/mingw64/include", + "-IC:/rtools40/mingw64/lib/gcc/x86_64-w64-mingw32/8.3.0/include-fixed", + "-IC:/rtools40/mingw64/x86_64-w64-mingw32/include" + ], + "pch": "RcppArmadillo", + "is_cpp": true, + "hash": "165330913916533285431653328560", + "compiler": "2977220399556355138", + "rversion": "4.1.2" +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/cpp-definition-cache b/.Rproj.user/864BD013/cpp-definition-cache new file mode 100644 index 0000000..55a113a --- /dev/null +++ b/.Rproj.user/864BD013/cpp-definition-cache @@ -0,0 +1,161 @@ +[ + { + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "file_last_write": 1653248994.0, + "definitions": [ + { + "usr": "c:@F@grad#$@N@Rcpp@S@Function_Impl>#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#d#", + "kind": 6, + "parent_name": "", + "name": "grad(Rcpp::Function, Rcpp::NumericVector, double)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 11, + "column": 21 + }, + { + "usr": "c:@F@binCoef#I#I#", + "kind": 6, + "parent_name": "", + "name": "binCoef(int, int)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 25, + "column": 5 + }, + { + "usr": "c:@F@gradN#$@N@Rcpp@S@Function_Impl>#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#d#I#", + "kind": 6, + "parent_name": "", + "name": "gradN(Rcpp::Function, Rcpp::NumericVector, double, int)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 32, + "column": 21 + }, + { + "usr": "c:@F@reltol#d#d#", + "kind": 6, + "parent_name": "", + "name": "reltol(double, double)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 74, + "column": 8 + }, + { + "usr": "c:@F@order_#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "order_(Rcpp::NumericVector)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 78, + "column": 15 + }, + { + "usr": "c:@F@mean_v#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#I#", + "kind": 6, + "parent_name": "", + "name": "mean_v(Rcpp::NumericVector, int)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 84, + "column": 8 + }, + { + "usr": "c:@F@crit_v#$@N@Rcpp@S@Matrix>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "crit_v(Rcpp::NumericMatrix)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 96, + "column": 8 + }, + { + "usr": "c:@F@steepSGD#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI19#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "steepSGD(Rcpp::NumericVector, Rcpp::List)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 109, + "column": 8 + }, + { + "usr": "c:@F@steepSTGD#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI19#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "steepSTGD(Rcpp::NumericVector, Rcpp::List)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 117, + "column": 8 + }, + { + "usr": "c:@F@steepLMM#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI19#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "steepLMM(Rcpp::NumericVector, Rcpp::List)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 130, + "column": 8 + }, + { + "usr": "c:@F@fibonacci#d#d#$@N@Rcpp@S@Vector>#VI19#@N@Rcpp@ST>1#T@PreserveStorage#I#", + "kind": 6, + "parent_name": "", + "name": "fibonacci(double, double, Rcpp::List, int)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 141, + "column": 8 + }, + { + "usr": "c:Algorithms.cpp@T@funcPtr", + "kind": 8, + "parent_name": "", + "name": "funcPtr", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 173, + "column": 18 + }, + { + "usr": "c:@F@putFunPtrInXPtr#$@N@std@N@__cxx11@S@basic_string>#C#$@N@std@S@char_traits>#C#$@N@std@S@allocator>#C#", + "kind": 6, + "parent_name": "", + "name": "putFunPtrInXPtr(std::string)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 175, + "column": 15 + }, + { + "usr": "c:@F@nelder_mead#$@N@std@N@__cxx11@S@basic_string>#C#$@N@std@S@char_traits>#C#$@N@std@S@allocator>#C#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI19#@N@Rcpp@ST>1#T@PreserveStorage#", + "kind": 6, + "parent_name": "", + "name": "nelder_mead(std::string, Rcpp::NumericVector, Rcpp::List)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 186, + "column": 12 + }, + { + "usr": "c:@F@SGD#$@N@Rcpp@S@Function_Impl>#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#S0_#d#I#d#b#", + "kind": 6, + "parent_name": "", + "name": "SGD(Rcpp::Function, Rcpp::NumericVector, Rcpp::Function, double, int, double, bool)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 305, + "column": 6 + }, + { + "usr": "c:@F@STGD#$@N@Rcpp@S@Function_Impl>#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#S0_#d#I#d#b#", + "kind": 6, + "parent_name": "", + "name": "STGD(Rcpp::Function, Rcpp::NumericVector, Rcpp::Function, double, int, double, bool)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 382, + "column": 6 + }, + { + "usr": "c:@F@LMM#$@N@Rcpp@S@Function_Impl>#@N@Rcpp@ST>1#T@PreserveStorage#$@N@Rcpp@S@Vector>#VI14#@N@Rcpp@ST>1#T@PreserveStorage#S0_#d#I#d#b#", + "kind": 6, + "parent_name": "", + "name": "LMM(Rcpp::Function, Rcpp::NumericVector, Rcpp::Function, double, int, double, bool)", + "file": "C:/Users/vrfra/Documents/GitHub/optimg/src/Algorithms.cpp", + "line": 463, + "column": 6 + } + ] + } +] \ No newline at end of file diff --git a/.Rproj.user/864BD013/pcs/debug-breakpoints.pper b/.Rproj.user/864BD013/pcs/debug-breakpoints.pper new file mode 100644 index 0000000..4893a8a --- /dev/null +++ b/.Rproj.user/864BD013/pcs/debug-breakpoints.pper @@ -0,0 +1,5 @@ +{ + "debugBreakpointsState": { + "breakpoints": [] + } +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/pcs/files-pane.pper b/.Rproj.user/864BD013/pcs/files-pane.pper new file mode 100644 index 0000000..67280ec --- /dev/null +++ b/.Rproj.user/864BD013/pcs/files-pane.pper @@ -0,0 +1,9 @@ +{ + "sortOrder": [ + { + "columnIndex": 2, + "ascending": true + } + ], + "path": "~/GitHub/optimg" +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/pcs/source-pane.pper b/.Rproj.user/864BD013/pcs/source-pane.pper new file mode 100644 index 0000000..a528f3b --- /dev/null +++ b/.Rproj.user/864BD013/pcs/source-pane.pper @@ -0,0 +1,3 @@ +{ + "activeTab": -1 +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/pcs/windowlayoutstate.pper b/.Rproj.user/864BD013/pcs/windowlayoutstate.pper new file mode 100644 index 0000000..544a777 --- /dev/null +++ b/.Rproj.user/864BD013/pcs/windowlayoutstate.pper @@ -0,0 +1,14 @@ +{ + "left": { + "splitterpos": 376, + "topwindowstate": "HIDE", + "panelheight": 904, + "windowheight": 942 + }, + "right": { + "splitterpos": 565, + "topwindowstate": "NORMAL", + "panelheight": 904, + "windowheight": 942 + } +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/pcs/workbench-pane.pper b/.Rproj.user/864BD013/pcs/workbench-pane.pper new file mode 100644 index 0000000..f398270 --- /dev/null +++ b/.Rproj.user/864BD013/pcs/workbench-pane.pper @@ -0,0 +1,5 @@ +{ + "TabSet1": 3, + "TabSet2": 3, + "TabZoom": {} +} \ No newline at end of file diff --git a/.Rproj.user/864BD013/persistent-state b/.Rproj.user/864BD013/persistent-state new file mode 100644 index 0000000..b304cff --- /dev/null +++ b/.Rproj.user/864BD013/persistent-state @@ -0,0 +1,8 @@ +build-last-errors="[]" +build-last-errors-base-dir="~/GitHub/optimg/" +build-last-outputs="[{\"type\":0,\"output\":\"==> Rcpp::compileAttributes()\\n\\n\"},{\"type\":1,\"output\":\"* Updated R/RcppExports.R\\r\\n\"},{\"type\":1,\"output\":\"\\n\"},{\"type\":0,\"output\":\"==> Rcmd.exe INSTALL --no-multiarch --with-keep.source optimg\\n\\n\"},{\"type\":1,\"output\":\"* installing to library 'C:/Users/vrfra/AppData/Local/R/win-library/4.2'\\r\\n\"},{\"type\":1,\"output\":\"* installing *source* package 'optimg' ...\\r\\n\"},{\"type\":1,\"output\":\"** using staged installation\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"make: Nothing to be done for 'all'.\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** libs\\r\\n\"},{\"type\":1,\"output\":\"installing to C:/Users/vrfra/AppData/Local/R/win-library/4.2/00LOCK-optimg/00new/optimg/libs/x64\\r\\n\"},{\"type\":1,\"output\":\"** R\\r\\n\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\r\\n\"},{\"type\":1,\"output\":\"*** installing help indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** building package indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\r\\n\"},{\"type\":1,\"output\":\"* DONE (optimg)\\r\\n\"},{\"type\":1,\"output\":\"\"}]" +compile_pdf_state="{\"tab_visible\":false,\"running\":false,\"target_file\":\"\",\"output\":\"\",\"errors\":[]}" +files.monitored-path="" +find-in-files-state="{\"handle\":\"\",\"input\":\"\",\"path\":\"\",\"regex\":false,\"ignoreCase\":false,\"results\":{\"file\":[],\"line\":[],\"lineValue\":[],\"matchOn\":[],\"matchOff\":[],\"replaceMatchOn\":[],\"replaceMatchOff\":[]},\"running\":false,\"replace\":false,\"preview\":false,\"gitFlag\":false,\"replacePattern\":\"\"}" +imageDirtyState="0" +saveActionState="0" diff --git a/.Rproj.user/864BD013/rmd-outputs b/.Rproj.user/864BD013/rmd-outputs new file mode 100644 index 0000000..3f2ff2d --- /dev/null +++ b/.Rproj.user/864BD013/rmd-outputs @@ -0,0 +1,5 @@ + + + + + diff --git a/.Rproj.user/864BD013/saved_source_markers b/.Rproj.user/864BD013/saved_source_markers new file mode 100644 index 0000000..2b1bef1 --- /dev/null +++ b/.Rproj.user/864BD013/saved_source_markers @@ -0,0 +1 @@ +{"active_set":"","sets":[]} \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 150169c..8e10ddf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,27 @@ -Package: optimg Type: Package -Title: General-purpose Gradient-based Optimization -Version: 0.1.0 -Year: 2021, October -Author: Vithor Rosa Franco +Package: optimg +Title: General-Purpose Gradient-Based Optimization +Version: 0.2.0 +Date: 2022-05-19 +Authors@R: c(person(given = "Vithor Rosa", family = "Franco", role = c("aut", "cre", "cph"), email = "vithorfranco@gmail.com"), person(given="Alexandre Augusto de Sá", family = "dos Santos", role = c("aut", "cph")), person(given = "Felipe Ofugi", family = "Hara", role = c("aut", "cph"))) +Author: Vithor Rosa Franco [aut, cre, cph], + Alexandre Augusto de Sá dos Santos [aut, cph], + Felipe Ofugi Hara [aut, cph] Maintainer: Vithor Rosa Franco -Description: This package is a general purpose tool for helping users to implement - gradient descent methods for function optimization. Currently, the - Steepest 2-Groups Gradient Descent and the Adaptive Moment Estimation - (Adam) are the methods implemented. Other methods will be implemented - in the future. -Imports: ucminf (>= 1.1-4) +Description: Provides general purpose tools for helping users to implement steepest + gradient descent methods for function optimization; for details see + Ruder (2016) . Currently, the Stepeest Gradient + Descent, the Steepest 2-Groups Gradient Descent and the Adaptive Moment + Estimation (Adam) are the methods implemented. Other methods will be + implemented in the future. License: GPL-3 -Encoding: UTF-8 +Depends: R(>= 4.0) +Imports: Rcpp (>= 1.0.8) +LinkingTo: Rcpp, + RcppArmadillo +SystemRequirements: C++11 LazyData: true URL: https://github.com/vthorrf/optimg +BugReports: https://github.com/vthorrf/optimg/issues +Encoding: UTF-8 +NeedsCompilation: yes \ No newline at end of file diff --git a/NAMESPACE b/NAMESPACE index 56cbbd9..83ac22d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,4 @@ +useDynLib(optimg, .registration=TRUE) +importFrom(Rcpp, evalCpp) #exportPattern("^[[:alpha:]]+") -export(optimg) -importFrom("ucminf", "ucminf") -importFrom("stats", "optim", "rnorm", "sd") -importFrom("utils", "setTxtProgressBar", "txtProgressBar") +export(optimg, gradient) \ No newline at end of file diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..4e4e1b4 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,23 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +grad <- function(Model, par, h = 1e-6) { + .Call(`_optimg_grad`, Model, par, h) +} + +gradN <- function(f, par, h = 1e-6, order = 1L) { + .Call(`_optimg_gradN`, f, par, h, order) +} + +SGD <- function(fn, startvalue, gr, h = 1e-6, maxit = 10L, tol = 1e-8, verbose = FALSE) { + .Call(`_optimg_SGD`, fn, startvalue, gr, h, maxit, tol, verbose) +} + +STGD <- function(fn, startvalue, gr, h = 1e-6, maxit = 10L, tol = 1e-8, verbose = FALSE) { + .Call(`_optimg_STGD`, fn, startvalue, gr, h, maxit, tol, verbose) +} + +LMM <- function(fn, startvalue, gr, h = 1e-6, maxit = 10L, tol = 1e-8, verbose = FALSE) { + .Call(`_optimg_LMM`, fn, startvalue, gr, h, maxit, tol, verbose) +} + diff --git a/R/algorithms.R b/R/algorithms.R deleted file mode 100644 index 44df8fe..0000000 --- a/R/algorithms.R +++ /dev/null @@ -1,287 +0,0 @@ -### Numerical gradient==== -grad <- function(mod, par, ..., Interval=1e-6) { - mat <- matrix(par, nrow=length(par), ncol=length(par)) - for (i in 1:ncol(mat)) { - mat[i,i] <- par[i] + Interval - } - df <- vector("numeric", length=length(par)) - f_x <- mod(par, ...) - for (i in 1:ncol(mat)) { - df[i] <- (mod(mat[,i], ...) - f_x) / Interval - } - df[which(!is.finite(df))] <- 0 - return(df) - return(df) -} -steep <- function(mod, par, ..., est, df) { - mod(est - {par[1]*df*{(abs(df) > sd(df)) * 1}} - - {par[2]*df*{(abs(df) <= sd(df)) * 1}}, ...) -} -steepA <- function(mod, par, ..., est, df) { - mod(est - {par*df}, ...) -} -Gammas <- function(mod, par, ..., est, df, vi, si, iter, epsilon) { - v <- {{par[1] * vi} + {{1 - par[1]} * df }}/{1 - {par[1] ^ iter}} - s <- {{par[2] * si} + {{1 - par[2]} * {df*df}}}/{1 - {par[2] ^ iter}} - newV <- est - {{par[3]*v}/{epsilon+sqrt(s)}} - return(mod(newV, ...)) -} -sHAR <- function(par, mod, ...) { - theta <- rnorm(length(par)) - d <- theta / sqrt(sum(theta*theta)) - u <- suppressWarnings(optim(.5, steepA, mod=mod, df=d, ..., est=par, - method="Nelder-Mead")$par) - prop <- par - {u * d} - return(list("prop"=prop, "u"=u, "d"=d)) -} -reltol <- function(x, tol) tol * (abs(x) + tol) - -### Adam==== -ADAM <- function(fn, startvalue, gr, ..., Interval=1e-6, - maxit=100, tol=1e-8, verbose=T) { - # Opening message - if(verbose==T) { - cat("Steepest Adam will run for ", maxit, " iterations at most.\n\n", sep="") - startTime = proc.time() - } - - # Initial settings - epsilon <- Interval - if(verbose==T) { - pb <- txtProgressBar(min=0, max=maxit, style=3) - } - par <- s <- v <- G <- matrix(NA, nrow=maxit, ncol=length(startvalue)) - GG <- matrix(NA, nrow=maxit, ncol=3) - f <- vector("numeric", length=maxit) - convergence = F - - # First step - G[1,] <- gr(fn, startvalue, ..., Interval=Interval) - alpha <- suppressWarnings(unlist(optim(0, fn=steepA, mod=fn, df=G[1,], - ..., est=startvalue, - method="Nelder-Mead")$par)) - GG[1,] <- c(.5, .5, alpha) - v[1,] <- G[1,]*GG[1,1] - s[1,] <- G[1,]*GG[1,2] - par[1,] <- startvalue - {alpha*G[1,]} - f[1] <- fn(par[1,], ...) - if(verbose==T) { - setTxtProgressBar(pb, 1) - } - - # Run ADAM algorithm - for(i in 2:maxit) { - G[i,] <- gr(fn, par[i-1,], ..., Interval=Interval) - GG[i,] <- suppressWarnings(unlist(optim(par=c(0,0,0), fn=Gammas, mod=fn, - df=G[i,], ..., est=par[i-1,], - vi=v[i-1,], si=s[i-1,], iter=i, - epsilon=epsilon)$par)) - if (sum(abs(GG[i,])) == 0) { - GG[i,] <- suppressWarnings(unlist(ucminf(par=c(0,0,0), fn=Gammas, mod=fn, - df=G[i,], ..., est=par[i-1,], - vi=v[i-1,], si=s[i-1,], iter=i, - epsilon=epsilon)$par)) - } - gammav <- GG[i,1] - gammas <- GG[i,2] - alpha <- GG[i,3] - v[i,] <- {{gammav * v[i-1,]} + {{1 - gammav} * G[i,] }}/{1 - {gammav^i}} - s[i,] <- {{gammas * s[i-1,]} + {{1 - gammas} * {G[i,]*G[i,]}}}/{1 - {gammas^i}} - par[i,] <- par[i-1,] - {{alpha*v[i,]}/{epsilon+sqrt(s[i,])}} - f[i] <- fn(par[i,], ...) - # Check convergence - if (reltol(f[i], tol) > abs(f[i] - f[i-1])) { - convergence = T - if (verbose==T) { - setTxtProgressBar(pb, maxit) - cat("\nConvergence achieved!") - } - break - } else { - if (verbose==T) { - setTxtProgressBar(pb, i) - } - } - } - if ({convergence==F} & {verbose==T}) { - cat("\nConvergence may not have been achieved!") - } - if (i < maxit) { - f <- f[-c({i+1}:maxit)] - par <- as.matrix(par[-c({i+1}:maxit),]) - G <- G[-c({i+1}:maxit),] - s <- s[-c({i+1}:maxit),] - v <- v[-c({i+1}:maxit),] - GG <- GG[-c({i+1}:maxit),] - } - if (verbose==T) { - close(pb) - cat("\n") - # Final messages - stopTime = proc.time() - elapsedTime = stopTime - startTime - cat("It took ",round(elapsedTime[3],2), - " secs for the run to finish. \n", sep="") - } - - # Return results - Results <- list("Cost"=f, "Estimates"=par, "Gradients"=G, - "DSG"=s, "Momentum"=v, "par"=par[which.min(f),], - "value"=f[which.min(f)], "steps"=GG) - return(Results) -} - -### Steepest 2-group Gradient Descent==== -STDM <- function(fn, startvalue, gr, ..., Interval=1e-6, maxit=100, tol=1e-8, - verbose=T) { - # Opening message - if(verbose==T) { - cat("Steepest 2-group Gradient Descent will run for ", - maxit, " iterations at most.\n\n", sep="") - startTime = proc.time() - } - - # Initial settings - if(verbose==T) { - pb <- txtProgressBar(min=0, max=maxit, style=3) - } - #par <- df <- array(dim = c(maxit,length(startvalue))) - par <- df <- matrix(NA, nrow=maxit, ncol=length(startvalue)) - f <- vector("numeric", length=maxit) - step <- array(dim = c(maxit,2)) - convergence = F - - # First step - df[1,] <- gr(fn, startvalue, ..., Interval=Interval) - step[1,] <- suppressWarnings(unlist(optim(c(0,0), fn=steep, mod=fn, - ..., df=df[1,], est=startvalue, - method="Nelder-Mead")$par)) - par[1,] <- startvalue - - {step[1,1]*df[1,]*{(abs(df[1,]) > sd(df[1,])) * 1}} - - {step[1,2]*df[1,]*{(abs(df[1,]) <= sd(df[1,])) * 1}} - f[1] <- fn(par[1,], ...) - if(verbose==T) { - setTxtProgressBar(pb, 1) - } - - # Start estimation - for(run in 2:maxit) { - # Calculate gradient and estimate step parameter - df[run,] <- gr(fn, par[run-1,], ..., Interval=Interval) - step[run,] <- suppressWarnings(unlist(optim(step[run-1,], fn=steep, mod=fn, - ..., df=df[run,], est=par[run-1,], - method="Nelder-Mead")$par)) - if (sum(abs(step[run,])) == 0) { - step[run,] <- suppressWarnings(unlist(ucminf(par=c(0,0), fn=steep, mod=fn, - ..., df=df[run,], - est=par[run-1,])$par)) - } - par[run,] <- par[run-1,] - - {step[run,1]*df[run,]*{(abs(df[run,]) > sd(df[run,])) * 1}} - - {step[run,2]*df[run,]*{(abs(df[run,]) <= sd(df[run,])) * 1}} - f[run] <- fn(par[run,], ...) - # Check convergence - if (reltol(f[run], tol) > abs(f[run] - f[run-1])) { - convergence = T - if(verbose==T) { - setTxtProgressBar(pb, maxit) - cat("\nConvergence achieved!") - } - break - } else { - if(verbose==T) { - setTxtProgressBar(pb, run) - } - } - } - if ({convergence==F} && {verbose==T}) { - cat("\nConvergence may not have been achieved!") - } - if (run < maxit) { - f <- f[-c({run+1}:maxit)] - par <- par[-c({run+1}:maxit),] - df <- df[-c({run+1}:maxit),] - step <- step[-c({run+1}:maxit),] - } - if(verbose==T) { - close(pb) - cat("\n") - # Final messages - stopTime = proc.time() - elapsedTime = stopTime - startTime - cat("It took ",round(elapsedTime[3],2)," secs for the run to finish. \n", - sep="") - } - - # Return results - Results <- list("Cost"=f, "Estimates"=par, "Gradients"=df, - "par"=par[which.min(f),], "value"=f[which.min(f)], - "steps"=step) - return(Results) -} - -### Steepest Hit-and-Run Gradient Descent==== -HRGD <- function(fn, startvalue, gr, ..., Interval=1e-6, maxit=100) { - # Opening message - cat("Hit-and-Run Gradient Descent will run for ", - maxit, " iterations at most.\n\n", sep="") - startTime = proc.time() - - # Initial settings - pb <- txtProgressBar(min=0, max=maxit, style=3) - par <- df <- array(dim = c(maxit,length(startvalue))) - f <- vector("numeric", length=maxit) - step <- vector("numeric", length=maxit) - convergence = F - - # First step - df[1,] <- gr(fn, startvalue, ..., Interval=Interval) - SS <- suppressWarnings(unlist(optim(c(0,0), fn=steep, mod=fn, - ..., df=df[1,], est=startvalue, - method="Nelder-Mead")$par)) - step[1] <- mean(SS) - par[1,] <- startvalue - - {SS[1]*df[1,]*{(abs(df[1,]) > sd(df[1,])) * 1}} - - {SS[2]*df[1,]*{(abs(df[1,]) <= sd(df[1,])) * 1}} - f[1] <- fn(par[1,], ...) - setTxtProgressBar(pb, 1) - - # Start estimation - for(run in 2:maxit) { - # Calculate gradient and estimate step parameter - temp <- sHAR(par=par[run-1,], mod=fn, ...) - df[run,] <- temp$d - step[run] <- temp$u - par[run,] <- temp$prop - f[run] <- fn(par[run,], ...) - # Check convergence - if ({step[run] - step[run-1]} == 0) { - convergence = T - setTxtProgressBar(pb, maxit) - cat("\nConvergence achieved!") - break - } else { setTxtProgressBar(pb, run) } - } - if (convergence == F) { - cat("\nConvergence may not have been achieved!") - } - if (run < maxit) { - f <- f[-c({run+1}:maxit)] - par <- par[-c({run+1}:maxit),] - df <- df[-c({run+1}:maxit),] - step <- step[-c({run+1}:maxit)] - } - close(pb) - cat("\n") - - # Final messages - stopTime = proc.time() - elapsedTime = stopTime - startTime - cat("It took ",round(elapsedTime[3],2)," secs for the run to finish. \n", sep="") - - # Return results - Results <- list("Cost"=f, "Estimates"=par, "Gradients"=df, - "par"=par[which.min(f),], "value"=f[which.min(f)], "steps"=step) - return(Results) -} - diff --git a/R/gradient.R b/R/gradient.R new file mode 100644 index 0000000..9dc828b --- /dev/null +++ b/R/gradient.R @@ -0,0 +1,5 @@ +# Wrapper +gradient <- function(par, fn, ..., Interval=1e-6, order=1) { + FN <- function(par) fn(par, ...) + return(gradN(f=FN, par=par, h=Interval, order=order)) +} diff --git a/R/optimg.R b/R/optimg.R index 0836b41..76e6862 100644 --- a/R/optimg.R +++ b/R/optimg.R @@ -1,33 +1,44 @@ # Wrapper -optimg <- function(par, fn, gr=NULL, ..., method=c("STGD","ADAM"), +optimg <- function(par, fn, gr=NULL, ..., method=c("SGD","STGD","LMM","ADAM"), Interval=1e-6, maxit=100, tol=1e-8, full=F, verbose=F) { ### Initial checks - if(is.null(gr)) { - gr <- grad - } - if(is.vector(method)) { + ## Method + if(is.vector(method) | is.null(method)) { method = method[1] } if({length(par)==1} & {method=="STGD"}) { - method="ADAM" + method="SGD" } - + ## Cost function + FN <- function(par) fn(par, ...) + ## Gradient function + if (!is.null(gr)) { + GR <- function(par) gr(par, ...) + } else if(is.null(gr)) { + GR <- function(par) grad(FN, par, Interval) + } else stop("Specification of gradient function is incorrect.") + ### Select method - if(method=="STGD"){ - Result <- STDM(fn, par, gr, ..., Interval=Interval, maxit=maxit, tol=tol, - verbose=verbose) + if(method=="SGD"){ + Result <- SGD(FN, par, GR, Interval, maxit, tol, verbose) + } else if(method=="STGD"){ + Result <- STGD(FN, par, GR, Interval, maxit, tol, verbose) + } else if(method=="LMM"){ + Result <- LMM(FN, par, GR, Interval, maxit, tol, verbose) } else if(method=="ADAM"){ - Result <- ADAM(fn, par, gr, ..., Interval=Interval, maxit=maxit, tol=tol, - verbose=verbose) - } else stop("Unkown method. Please, see optimg's documentation.") - + Result <- ADAM(FN, par, GR, Interval, maxit, tol, verbose) + } else stop("Unkown method. Please check optimg's documentation.") + ### Return results + par <- as.matrix(Result$Estimates)[Result$MaxInt,] + value <- Result$Cost[Result$MaxInt] + counts <- Result$MaxInt + convergence <- Result$convergence if(full==T) { - return(Result) + return(list("par"=par, "value"=value, "counts"=counts, + "convergence"=convergence, "Full"=Result)) } else if(full==F) { - Results <- list("par"=Result$par, "value"=Result$value, - "counts"=length(Result$steps), - "convergence"={length(Result$steps) == maxit} * 1) + Results <- list("par"=par, "value"=value, "counts"=counts, "convergence"=convergence) return(Results) } else stop("Incorrect value for argument 'full'.") } diff --git a/README.md b/README.md index 210f94c..ca33c68 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ optimg: General-purpose Gradient-based Optimization (version 0.1.2) CRAN checks

-This package is a general purpose tool for helping users to implement gradient descent methods for function optimization. Currently, the Steepest 2-Groups Gradient Descent and the Adaptive Moment Estimation (Adam) are the methods implemented. Other methods will be implemented in the future. +This package is a general purpose tool for helping users to implement gradient descent methods for function optimization. Currently, the Steepest Gradient Descent, the Steepest 2-Groups Gradient Descent, the Levenberg–Marquardt algorithm, and the Adaptive Moment Estimation (Adam) are the methods implemented. Other methods will be implemented in the future. This package should be considered experimental in this point of development. diff --git a/man/gradient.Rd b/man/gradient.Rd new file mode 100644 index 0000000..1628f94 --- /dev/null +++ b/man/gradient.Rd @@ -0,0 +1,43 @@ +\name{gradient} +\alias{gradient} +\title{Numerical Gradient of a Function} +\description{ +Calculate the gradient of a function by finite difference coefficient. +} +\usage{ +gradient(par, fn, ..., Interval=1e-6, order=1) +} +\arguments{ + \item{par}{Initial values for the parameters to be optimized over.} + \item{fn}{A function to be minimized. It should return a scalar result.} + \item{...}{Further arguments to be passed to fn.} + \item{Interval}{The epsilon difference to be used by gr.} + \item{order}{An integer expressing the order of the derivative.} +} +\details{ +This function calculates a numerical approximation of the kth-order derivative of the supplied function at the points par. The arguments in ... are passed to fn, and it is assumed that fn returns a single value. +} +\value{ +A real scalar or vector of the approximated gradient(s). +} +\examples{ +### Fit a simple linear regression with RMSE as cost function +require(optimg) + +# Predictor +x <- seq(-3,3,len=100) +# Criterion +y <- rnorm(100, 2 + {1.2*x}, 1) + +# RMSE cost function +fn <- function(par, Y, X) { + mu <- par[1] + {par[2] * X} + rmse <- sqrt(mean({Y-mu}^2)) + return(rmse) +} + +# Return gradients +gradient(c(0,0), fn, Y=y, X=x, order=1) +gradient(c(0,0), fn, Y=y, X=x, order=2) +gradient(c(0,0), fn, Y=y, X=x, order=3) +} diff --git a/man/optimg.Rd b/man/optimg.Rd index 0feed8f..f465fde 100644 --- a/man/optimg.Rd +++ b/man/optimg.Rd @@ -1,17 +1,17 @@ \name{optimg} \alias{optimg} -\title{General-purpose Gradient-based Optimization} +\title{General-Purpose Gradient-Based Optimization} \description{ General-purpose optimization based on gradient descent algorithms. It is greatly inspired on the stats::optim function, aiming at increased usability and adaptability for optim's users. } \usage{ -optimg(par, fn, gr=NULL, ..., method=c("STGD","ADAM"), +optimg(par, fn, gr=NULL, ..., method=c("SGD","STGD","LMM","ADAM"), Interval=1e-6, maxit=100, tol=1e-8, full=F, verbose=F) } \arguments{ \item{par}{Initial values for the parameters to be optimized over.} \item{fn}{A function to be minimized. It should return a scalar result.} - \item{gr}{A function to return the gradient. Should include "Interval" as an argument, even if not used. If it is NULL, a finite-difference approximation will be used.} + \item{gr}{A function to return the gradient. If a function is not provided, a finite-difference approximation will be used.} \item{...}{Further arguments to be passed to fn and gr.} \item{method}{The method to be used. See ‘Details’.} \item{Interval}{The epsilon difference to be used by gr.} @@ -23,7 +23,11 @@ optimg(par, fn, gr=NULL, ..., method=c("STGD","ADAM"), \details{ As with the optim function, ... must be matched exactly. Also as with optim, optimg performs minimization. All the methods implemented are the "steepest" variation of the original methods. This means that tuning parameters are optimized at each step automatically. This makes the algorithms slower, but also more effective, especially for complex models. -The default method (unless the length of par is equal to 1; in which case, the default is "ADAM") is an implementation of the Steepest 2-Group Gradient Descent algorithm. This algorithm is a variation of the Steepest Gradient Descent method which calculates different step sizes for two groups of gradients: those within a standard deviation (below or above), and those beyond a standard deviation (below or above). +The default method is the traditional Steepest Gradient Descent ("SGD") algorithm. + +Method "STGD" is a variation of the Steepest Gradient Descent method which optimizes different step sizes for two groups of gradients: those within a standard deviation (below or above the mean), and those beyond a standard deviation (below or above the mean). If the specified function has only one parameter, SGD is used instead. + +Method "LMM" is the Levenberg-Marquardt method. This method computes an approximate regularized truncated Hessian at each step. The LMM is a robust version of the Gaussian-Newton method, which means that it can better handle far off initial values. Method "ADAM" is the Adaptive Moment Estimation method. This method computes adaptive learning rates for each parameter. Adam stores both exponentially decaying average of past squared gradients, as well as measures of momentum. } @@ -45,17 +49,19 @@ x <- seq(-3,3,len=100) y <- rnorm(100, 2 + {1.2*x}, 1) # RMSE cost function -fn <- function(par, X) { +fn <- function(par, Y, X) { mu <- par[1] + {par[2] * X} - rmse <- sqrt(mean({y-mu}^2)) + rmse <- sqrt(mean({Y-mu}^2)) return(rmse) } # Compare optimization methods -optim(c(0,0),fn,X=x,method="Nelder-Mead") -optim(c(0,0),fn,X=x,method="BFGS") -optim(c(0,0),fn,X=x,method="CG") -optim(c(0,0),fn,X=x,method="L-BFGS-B") -optimg(c(0,0),fn,X=x,method="ADAM") -optimg(c(0,0),fn,X=x,method="STGD") +optim(c(0,0), fn, Y=y, X=x, method="Nelder-Mead") +optim(c(0,0), fn, Y=y, X=x, method="BFGS") +optim(c(0,0), fn, Y=y, X=x, method="CG") +optim(c(0,0), fn, Y=y, X=x, method="L-BFGS-B") +optimg(c(0,0), fn, Y=y, X=x, method="SGD") +optimg(c(0,0), fn, Y=y, X=x, method="STGD") +optimg(c(0,0), fn, Y=y, X=x, method="LMM") +optimg(c(0,0), fn, Y=y, X=x, method="ADAM") } diff --git a/optimg.Rproj b/optimg.Rproj index 497f8bf..21a4da0 100644 --- a/optimg.Rproj +++ b/optimg.Rproj @@ -12,9 +12,6 @@ Encoding: UTF-8 RnwWeave: Sweave LaTeX: pdfLaTeX -AutoAppendNewline: Yes -StripTrailingWhitespace: Yes - BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/src/Algorithms.cpp b/src/Algorithms.cpp new file mode 100644 index 0000000..7592590 --- /dev/null +++ b/src/Algorithms.cpp @@ -0,0 +1,545 @@ +// [[Rcpp::depends(RcppArmadillo)]] +#include +// [[Rcpp::plugins(cpp11)]] +#include + +using namespace Rcpp; +using namespace arma; + +// ====================================Supporting functions==================================== +// [[Rcpp::export]] +Rcpp::NumericVector grad(Function Model, NumericVector par, double h=1e-6) { + NumericMatrix mat(par.length(), par.length()); + NumericVector out(par.length()); + double f_x = as(Model(par)); + + for (int i = 0; i < mat.ncol(); i++) { + mat(i, _) = par; + mat(i, i) = mat(i, i) + h; + out[i] = (as(Model(mat(i,_))) - f_x) / h; + } + + return out; +} + +int binCoef(int n, int k) { + if (k == 0 || k == n) + return 1; + return binCoef(n - 1, k - 1) + binCoef(n - 1, k); +} + +// [[Rcpp::export]] +Rcpp::NumericVector gradN(Function f, NumericVector par, double h=1e-6, int order=1) { + int m = par.length(); + NumericMatrix df(m,order); + IntegerVector posit(m); + for(int i = 0; i < m; ++i) { + posit[i] = i; + } + for (int i = 0; i < m; ++i) { + for(int j = 0; j < order; ++j) { + NumericVector temp(m); + for(int k = 0; k < m; ++k) { + double t0; + if (posit[k] == i) { + t0 = 1.0; + } else { + t0 = 0.0; + } + temp[k] = par[k] + (t0 * double(j+1) * h); + } + df(i,j) = pow(-1.0, double(order-j+1)) * double(binCoef(order,j+1)) * as(f(temp)); + if(!arma::is_finite(df(i,j))) { + df(i,j) = 0; + } + } + } + + NumericMatrix Final(m,order+1); + double fix = as(f(par)) * pow(-1, order % 2); + for(int i = 0; i < m; ++i) { + Final(i,0) = fix; + } + for(int j = 1; j < (order+1); ++j) { + Final(_,j) = df(_,(j-1)); + } + NumericVector out(m); + for(int i = 0; i < m; ++i) { + out[i] = sum(Final(i,_))/h; + } + + return out; +} + +double reltol(double cost, double tol=1e-8) { + return tol * (abs(cost) + tol); +} + +IntegerVector order_(NumericVector x) { + // Picking up order() function from base package + Function f("order"); + return as(f(x))-1; +} + +double mean_v(NumericVector vec, int ele) { + int vs = vec.size(); + double res = 0; + for (int i = 0; i < vs; i++) { + if (i == ele) { + continue; + } + res += vec[i]; + } + return res / (vs-1); +} + +double crit_v(NumericMatrix Obj) { + int ms = Obj.cols(); + NumericVector x1(ms); + NumericMatrix matx(ms, ms); + for (int iter = 0; iter < ms; iter++) { + for (int j = 0; j < ms; j++) { + matx(iter, j) = abs(Obj(iter + 1, j) - Obj(iter, j)); + } + x1[iter] = max(matx(iter, _)); + } + return max(x1); +} + +double steepSGD(NumericVector par, List Data) { + Rcpp::NumericVector der = Data["Gradients"]; + Rcpp::NumericVector init = Data["Estimates"]; + Rcpp::Function FUN = Data["FUN"]; + + return as(FUN(init - (par * der))); +} + +double steepSTGD(NumericVector par, List Data) { + Rcpp::NumericVector der = Data["Gradients"]; + Rcpp::NumericVector init = Data["Estimates"]; + Rcpp::Function FUN = Data["FUN"]; + Rcpp::NumericVector New(init.length()); + double SDD = sd(der); + for (int iter = 0; iter < New.length(); iter++) { + New[iter] = abs(der[iter]) > SDD ? init[iter] - (par[0] * der[iter]) : init[iter] - (par[1] * der[iter]); + } + + return as(FUN(New)); +} + +double steepLMM(NumericVector par, List Data) { + arma::vec der = as(wrap(Data["Gradients"])); + Rcpp::NumericVector init = Data["Estimates"]; + Rcpp::Function FUN = Data["FUN"]; + Rcpp::NumericMatrix I = NumericMatrix::diag(init.length(), exp(par[1])); + arma::mat SOG = (der * der.t()) + as(wrap(I)); + Rcpp::NumericVector LM = as(wrap(arma::inv(SOG) * der)); + + return as(FUN(init - (par[0] * LM))); +} + +double fibonacci(double lower, double upper, List Data, int n=100) { + double epsilon = .01; + double a = lower; + double b = upper; + double phi = (1.0 + sqrt(5.0)) / 2.0; + double S = (1.0 - sqrt(5.0)) / phi; + double rho = 1.0 / (phi * (1 - pow(S, n + 1)) / (1 - pow(S, n))); + double d = (rho * b) + ((1 - rho) * a); + double yd = steepSGD(d, Data); + double c, yc; + for (int i = 0; i < n-1; i++) { + if(a == b) { + break; + } + if (i == (n - 1)) { + c = (epsilon * a) + ((1 - epsilon) * d); + } + else { + c = (rho * a) + ((1 - rho) * b); + } + yc = steepSGD(c, Data); + if (yc < yd) { + b = d, d = c, yd = yc; + } + else { + a = b, b = c; + } + rho = 1.0 / (phi * (1 - pow(S, n - i + 1)) / (1 - pow(S, n - i)));; + } + return (a + b) / 2.0; +} + +typedef double (*funcPtr)(NumericVector par, List Data); + +XPtr putFunPtrInXPtr(std::string fstr) { + if (fstr == "steepSGD") + return XPtr(new funcPtr(&steepSGD)); + else if (fstr == "steepSTGD") + return(XPtr(new funcPtr(&steepSTGD))); + else if (fstr == "steepLMM") + return(XPtr(new funcPtr(&steepLMM))); + else + return XPtr(R_NilValue); // runtime error as NULL no XPtr +} + +Rcpp::List nelder_mead(std::string funname, NumericVector init, List Data) { + XPtr xpfun = putFunPtrInXPtr(funname); + funcPtr fn = *xpfun; + double epsilon = 1e-10; + int d = init.length(); + int d1 = d + 1; + double alpha, beta, gamma, sigma; + int maxfeval = 50 * pow(d, 2); + if (d > 1) { + alpha = 1.0, beta = 1.0 + 2.0 / d, gamma = .75 - (.5 / d), sigma = 1.0 - (1.0/d); + } + else { + alpha = 1.0, beta = 2.0, gamma = .50, sigma=.50; + } + double xi = min(NumericVector::create(max(NumericVector::create(max(abs(init)), 1.0)), 10.0)); + + // Initial settings + NumericMatrix S(d1, d); + for (int iter = 0; iter < d1; iter++) { + for (int j = 0; j < d; j++) { + S(iter, j) = iter == j ? 1.0 : iter == d ? (1.0 - sqrt(d1))/d : 0.0; + } + } + NumericMatrix X(d1, d); + NumericVector y_arr(d1); + for (int iter = 0; iter < d1; iter++) { + X(iter,_) = init + xi * S(iter,_); + y_arr[iter] = fn(X(iter,_), Data); + } + // Ordering + IntegerVector p = order_(y_arr); + NumericMatrix SI(d1, d); + NumericVector y_I(d1); + for (int iter = 0; iter < d1; iter++) { + SI(iter,_) = X(p[iter],_); + y_I[iter] = y_arr[p[iter]]; + } + X = SI; + y_arr = y_I; + int ct = d1; + double delta = crit_v(X); + + while (delta >= (xi * epsilon)) { + if (ct > maxfeval) { + break; + } + // Simplex + NumericVector xl = X(0,_); + double yl = y_arr[0]; + NumericVector xh = X(d,_); + double yh = y_arr[d]; + NumericVector xs = X(d-1, _); + double ys = y_arr[d-1]; + // Centroid + NumericVector xm(d); + for (int iter = 0; iter < d; iter++) { + xm[iter] = mean_v(X(_,iter), d); + } + // Reflection + NumericVector xr = xm + alpha * (xm - xh); + double yr = fn(xr, Data); + ct += 1; + + // Expansion + if (yr < yl) { + NumericVector xe = xm + beta * (xr - xm); + double ye = fn(xe, Data); + X(d, _) = ye < yr ? xe : xr; + y_arr[d] = ye < yr ? ye : yr; + ct += 1; + } + else if (yr >= ys) { + if (!(yr >= yh)) { + xh = xr; + yh = yr; + X(d, _) = xr; + y_arr[d] = yr; + } + // Contraction + NumericVector xc = xm + gamma * (xh - xm); + double yc = fn(xc, Data); + ct += 1; + // Shrinkage + if (yc > yh) { + for (int iter = 1; iter < d1; iter++) { + X(iter, _) = xl + sigma * (X(iter, _) - xl); + y_arr[iter] = fn(X(iter, _), Data); + } + ct += d; + } + else { + X(d,_) = xc; + y_arr[d] = yc; + } + } + else { + X(d, _) = xr; + y_arr[d] = yr; + } + p = order_(y_arr); + NumericMatrix SI(d1, d); + NumericVector y_I(d1); + for (int iter = 0; iter < d1; iter++) { + SI(iter, _) = X(p[iter], _); + y_I[iter] = y_arr[p[iter]]; + } + X = SI; + y_arr = y_I; + delta = crit_v(X); + } + + return wrap(Rcpp::List::create(Rcpp::Named("par") = as(wrap(X(0,_))), + Rcpp::Named("value") = y_arr[0], + Rcpp::Named("nfeval") = ct, + Rcpp::Named("convergence") = delta > epsilon)); +} + +// ====================================Gradient Descent Algorithms==================================== +// [[Rcpp::export]] +SEXP SGD(Function fn, NumericVector startvalue, Function gr, double h=1e-6, + int maxit=10, double tol=1e-8, bool verbose=false) { + // Opening message + if (verbose == true) { + Rcpp::Rcout << "Steepest Gradient Descent will run for " << maxit << " iterations at most." << std::endl; + } + + // Initial settings + if(h > tol) { + tol = h; + } + NumericMatrix par(maxit, startvalue.length()); + NumericMatrix df(maxit, startvalue.length()); + NumericVector f(maxit); + NumericVector step(maxit); + bool convergence = false; + int final = 0; + + // First step + df(0, _) = as(gr(startvalue)); + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = startvalue, + Rcpp::Named("Gradients") = df(0, _), + Rcpp::Named("FUN") = fn); + step[0] = fibonacci(-1.0, 1.0, DataSGD, 1000); + NumericVector sss(1); + sss[0] = step[0]; + par(0,_) = startvalue - (step[0] * df(0,_)); + f[0] = as(fn(par(0, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << 1 << ". Cost: " << f[0] << std::endl; + } + + // Next steps + for (int run = 1; run < maxit; run++) { + // Gradients + df(run, _) = as(gr(par(run-1, _))); + // Steepest step + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = par(run-1,_), + Rcpp::Named("Gradients") = df(run, _), + Rcpp::Named("FUN") = fn); + //step[run] = as(wrap(nelder_mead("steepSGD", sss, DataSGD)["par"])); + step[run] = fibonacci(-1.0, 1.0, DataSGD, 1000); + // New parameter + par(run, _) = par(run-1,_) - (step[run] * df(run, _)); + f[run] = as(fn(par(run, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << run + 1 << ". Cost: " << f[run] << std::endl; + } + final += 1; + // Check convergence + if ( reltol(f[run], tol) > abs(f[run] - f[run-1]) ) { + convergence = true; + break; + } + } + if (verbose == true) { + if (convergence == true) { + Rcpp::Rcout << "Convergence achieved!" << std::endl; + } + else { + Rcpp::Rcout << "Convergence may not have been achieved!" << std::endl; + } + } + f = f[Range(0,final)]; + par = par(Range(0,final),_); + df = df(Range(0,final),_); + step = step[Range(0,final)]; + // Final Result + return wrap(Rcpp::List::create(Rcpp::Named("Cost") = f, + Rcpp::Named("Estimates") = par, + Rcpp::Named("Gradients") = df, + Rcpp::Named("steps") = step, + Rcpp::Named("MaxInt") = final + 1, + Rcpp::Named("convergence") = convergence)); +} + +// [[Rcpp::export]] +SEXP STGD(Function fn, NumericVector startvalue, Function gr, double h=1e-6, + int maxit=10, double tol=1e-8, bool verbose=false) { + // Opening message + if (verbose == true) { + Rcpp::Rcout << "Steepest 2-group Gradient Descent will run for " << maxit << + " iterations at most." << std::endl; + } + + // Initial settings + if(h > tol) { + tol = h; + } + NumericMatrix par(maxit, startvalue.length()); + NumericMatrix df(maxit, startvalue.length()); + NumericVector f(maxit); + NumericMatrix step(maxit, 2); + NumericVector sss = { .001,.001 }; + bool convergence = false; + int final = 0; + + // First step + df(0, _) = as(gr(startvalue)); + // Steepest step + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = startvalue, + Rcpp::Named("Gradients") = df(0, _), + Rcpp::Named("FUN") = fn); + step(0, _) = as(wrap(nelder_mead("steepSTGD", sss, DataSGD)["par"])); + for (int iter = 0; iter < startvalue.length(); iter++) { + par(0, iter) = abs(df(0, iter)) > sd(df(0, _)) ? startvalue[iter] - (step(0, 0) * df(0, iter)) : startvalue[iter] - (step(0, 1) * df(0, iter)); + } + f[0] = as(fn(par(0, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << 1 << ". Cost: " << f[0] << std::endl; + } + + // Next steps + for (int run = 1; run < maxit; run++) { + // Gradients + df(run, _) = as(gr(par(run - 1, _))); + // Steepest steps + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = par(run-1,_), + Rcpp::Named("Gradients") = df(run, _), + Rcpp::Named("FUN") = fn); + step(run, _) = as(wrap(nelder_mead("steepSTGD", sss, DataSGD)["par"])); + // New parameter + for (int iter = 0; iter < startvalue.length(); iter++) { + par(run, iter) = abs(df(run, iter)) > sd(df(run, _)) ? par(run - 1, iter) - (step(run, 0) * df(run, iter)) : par(run - 1, iter) - (step(run, 1) * df(run, iter)); + } + f[run] = as(fn(par(run, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << run + 1 << ". Cost: " << f[run] << std::endl; + } + final += 1; + // Check convergence + if ( reltol(f[run], tol) > abs(f[run] - f[run-1]) ) { + convergence = true; + break; + } + } + if (verbose == true) { + if (convergence == true) { + Rcpp::Rcout << "Convergence achieved!" << std::endl; + } + else { + Rcpp::Rcout << "Convergence may not have been achieved!" << std::endl; + } + } + f = f[Range(0,final)]; + par = par(Range(0,final),_); + df = df(Range(0,final),_); + step = step(Range(0,final),_); + // Final Result + return wrap(Rcpp::List::create(Rcpp::Named("Cost") = f, + Rcpp::Named("Estimates") = par, + Rcpp::Named("Gradients") = df, + Rcpp::Named("steps") = step, + Rcpp::Named("MaxInt") = final+1, + Rcpp::Named("convergence") = convergence)); +} + +// [[Rcpp::export]] +SEXP LMM(Function fn, NumericVector startvalue, Function gr, double h = 1e-6, + int maxit = 10, double tol = 1e-8, bool verbose = false) { + // Opening message + if (verbose == true) { + Rcpp::Rcout << "Levenberg-Marquardt method will run for " << maxit << " iterations at most." << std::endl; + } + + // Initial settings + if (h > tol) { + tol = h; + } + NumericMatrix par(maxit, startvalue.length()); + NumericMatrix df(maxit, startvalue.length()); + NumericVector f(maxit); + NumericMatrix step(maxit, 2); + NumericVector sss = { .001,.001 }; + bool convergence = false; + int final = 0; + + // First step + df(0, _) = as(gr(startvalue)); + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = startvalue, + Rcpp::Named("Gradients") = df(0, _), + Rcpp::Named("FUN") = fn); + //sss = fibonacci(-1.0, 1.0, DataSGD, 1000); + //par(0, _) = startvalue - (sss * df(0, _)); + step(0, _) = as(wrap(nelder_mead("steepLMM", sss, DataSGD)["par"])); + Rcpp::NumericMatrix I = NumericMatrix::diag(startvalue.length(), exp(step(0, 1))); + arma::mat SOG = (as(wrap(df(0, _))) * as(wrap(df(0, _))).t()) + as(wrap(I)); + Rcpp::NumericVector LM = as(wrap(arma::inv(SOG) * as(wrap(df(0, _))))); + par(0, _) = startvalue - (step(0, 0) * LM); + f[0] = as(fn(par(0, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << 1 << ". Cost: " << f[0] << std::endl; + } + + // Next steps + for (int run = 1; run < maxit; run++) { + // Gradients + df(run, _) = as(gr(par(run - 1, _))); + // Steepest step + List DataSGD = Rcpp::List::create(Rcpp::Named("Estimates") = par(run-1,_), + Rcpp::Named("Gradients") = df(run, _), + Rcpp::Named("FUN") = fn); + step(run, _) = as(wrap(nelder_mead("steepLMM", sss, DataSGD)["par"])); + // New parameter + Rcpp::NumericMatrix I = NumericMatrix::diag(startvalue.length(), exp(step(run,1))); + arma::mat SOG = (as(wrap(df(run, _))) * as(wrap(df(run, _))).t()) + as(wrap(I)); + Rcpp::NumericVector LM = as(wrap(arma::inv(SOG) * as(wrap(df(run, _))))); + par(run, _) = par(run-1, _) - (step(run,0) * LM); + f[run] = as(fn(par(run, _))); + if (verbose == true) { + Rcpp::Rcout << "Iteration: " << run + 1 << ". Cost: " << f[run] << std::endl; + } + final += 1; + // Check convergence + if (reltol(f[run], tol) > abs(f[run] - f[run - 1])) { + convergence = true; + break; + } + } + if (verbose == true) { + if (convergence == true) { + Rcpp::Rcout << "Convergence achieved!" << std::endl; + } + else { + Rcpp::Rcout << "Convergence may not have been achieved!" << std::endl; + } + } + f = f[Range(0, final)]; + par = par(Range(0, final), _); + df = df(Range(0, final), _); + step = step(Range(0, final), _); + // Final Result + return wrap(Rcpp::List::create(Rcpp::Named("Cost") = f, + Rcpp::Named("Estimates") = par, + Rcpp::Named("Gradients") = df, + Rcpp::Named("steps") = step, + Rcpp::Named("MaxInt") = final+1, + Rcpp::Named("convergence") = convergence)); +} + +// ====================================THE END==================================== \ No newline at end of file diff --git a/src/Algorithms.o b/src/Algorithms.o new file mode 100644 index 0000000..ab42bed Binary files /dev/null and b/src/Algorithms.o differ diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..bec5291 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,6 @@ +PKG_CPPFLAGS = -I../inst/include \ + -DBOOST_NO_LONG_LONG \ + -DBOOST_NO_AUTO_PTR \ + -DBOOST_BIND_GLOBAL_PLACEHOLDERS \ + -DRCPP_NO_RTTI \ + -DRCPP_USE_UNWIND_PROTECT \ No newline at end of file diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..bec5291 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,6 @@ +PKG_CPPFLAGS = -I../inst/include \ + -DBOOST_NO_LONG_LONG \ + -DBOOST_NO_AUTO_PTR \ + -DBOOST_BIND_GLOBAL_PLACEHOLDERS \ + -DRCPP_NO_RTTI \ + -DRCPP_USE_UNWIND_PROTECT \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..9d58c6f --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,105 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// grad +Rcpp::NumericVector grad(Function Model, NumericVector par, double h); +RcppExport SEXP _optimg_grad(SEXP ModelSEXP, SEXP parSEXP, SEXP hSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Function >::type Model(ModelSEXP); + Rcpp::traits::input_parameter< NumericVector >::type par(parSEXP); + Rcpp::traits::input_parameter< double >::type h(hSEXP); + rcpp_result_gen = Rcpp::wrap(grad(Model, par, h)); + return rcpp_result_gen; +END_RCPP +} +// gradN +Rcpp::NumericVector gradN(Function f, NumericVector par, double h, int order); +RcppExport SEXP _optimg_gradN(SEXP fSEXP, SEXP parSEXP, SEXP hSEXP, SEXP orderSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Function >::type f(fSEXP); + Rcpp::traits::input_parameter< NumericVector >::type par(parSEXP); + Rcpp::traits::input_parameter< double >::type h(hSEXP); + Rcpp::traits::input_parameter< int >::type order(orderSEXP); + rcpp_result_gen = Rcpp::wrap(gradN(f, par, h, order)); + return rcpp_result_gen; +END_RCPP +} +// SGD +SEXP SGD(Function fn, NumericVector startvalue, Function gr, double h, int maxit, double tol, bool verbose); +RcppExport SEXP _optimg_SGD(SEXP fnSEXP, SEXP startvalueSEXP, SEXP grSEXP, SEXP hSEXP, SEXP maxitSEXP, SEXP tolSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Function >::type fn(fnSEXP); + Rcpp::traits::input_parameter< NumericVector >::type startvalue(startvalueSEXP); + Rcpp::traits::input_parameter< Function >::type gr(grSEXP); + Rcpp::traits::input_parameter< double >::type h(hSEXP); + Rcpp::traits::input_parameter< int >::type maxit(maxitSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(SGD(fn, startvalue, gr, h, maxit, tol, verbose)); + return rcpp_result_gen; +END_RCPP +} +// STGD +SEXP STGD(Function fn, NumericVector startvalue, Function gr, double h, int maxit, double tol, bool verbose); +RcppExport SEXP _optimg_STGD(SEXP fnSEXP, SEXP startvalueSEXP, SEXP grSEXP, SEXP hSEXP, SEXP maxitSEXP, SEXP tolSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Function >::type fn(fnSEXP); + Rcpp::traits::input_parameter< NumericVector >::type startvalue(startvalueSEXP); + Rcpp::traits::input_parameter< Function >::type gr(grSEXP); + Rcpp::traits::input_parameter< double >::type h(hSEXP); + Rcpp::traits::input_parameter< int >::type maxit(maxitSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(STGD(fn, startvalue, gr, h, maxit, tol, verbose)); + return rcpp_result_gen; +END_RCPP +} +// LMM +SEXP LMM(Function fn, NumericVector startvalue, Function gr, double h, int maxit, double tol, bool verbose); +RcppExport SEXP _optimg_LMM(SEXP fnSEXP, SEXP startvalueSEXP, SEXP grSEXP, SEXP hSEXP, SEXP maxitSEXP, SEXP tolSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Function >::type fn(fnSEXP); + Rcpp::traits::input_parameter< NumericVector >::type startvalue(startvalueSEXP); + Rcpp::traits::input_parameter< Function >::type gr(grSEXP); + Rcpp::traits::input_parameter< double >::type h(hSEXP); + Rcpp::traits::input_parameter< int >::type maxit(maxitSEXP); + Rcpp::traits::input_parameter< double >::type tol(tolSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(LMM(fn, startvalue, gr, h, maxit, tol, verbose)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_optimg_grad", (DL_FUNC) &_optimg_grad, 3}, + {"_optimg_gradN", (DL_FUNC) &_optimg_gradN, 4}, + {"_optimg_SGD", (DL_FUNC) &_optimg_SGD, 7}, + {"_optimg_STGD", (DL_FUNC) &_optimg_STGD, 7}, + {"_optimg_LMM", (DL_FUNC) &_optimg_LMM, 7}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_optimg(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/RcppExports.o b/src/RcppExports.o new file mode 100644 index 0000000..07f3c1b Binary files /dev/null and b/src/RcppExports.o differ diff --git a/src/optimg.dll b/src/optimg.dll new file mode 100644 index 0000000..53b1812 Binary files /dev/null and b/src/optimg.dll differ