[Rcpp-commits] r2642 - in pkg/RcppDE/inst: . doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 2 05:59:12 CET 2010
Author: edd
Date: 2010-12-02 05:59:12 +0100 (Thu, 02 Dec 2010)
New Revision: 2642
Added:
pkg/RcppDE/inst/doc/
pkg/RcppDE/inst/doc/RcppDE.bib
pkg/RcppDE/inst/doc/RcppDE.tex
pkg/RcppDE/inst/doc/jss.bst
pkg/RcppDE/inst/doc/jss.cls
pkg/RcppDE/inst/doc/jsslogo.jpg
Log:
beginnings of a vignette for RcppDE
Added: pkg/RcppDE/inst/doc/RcppDE.bib
===================================================================
--- pkg/RcppDE/inst/doc/RcppDE.bib (rev 0)
+++ pkg/RcppDE/inst/doc/RcppDE.bib 2010-12-02 04:59:12 UTC (rev 2642)
@@ -0,0 +1,102 @@
+ at String{CRAN = "http://cran.r-project.org/" }
+ at String{manuals = CRAN # "doc/manuals/" }
+ at String{RCoreTeam = "{R Development Core Team}" }
+ at String{RFoundation = "R Foundation for Statistical Computing" }
+
+ at manual{R:exts,
+ author = RCoreTeam,
+ organization = RFoundation,
+ address = {Vienna, Austria},
+ year = 2010,
+ isbn = {3-900051-11-9},
+ title = "Writing R Extensions",
+ url = manuals # "R-exts.html"
+}
+
+ at Manual{R:Main,
+ title = {R: A Language and Environment for Statistical
+ Computing},
+ author = RCoreTeam,
+ organization = RFoundation,
+ address = {Vienna, Austria},
+ year = 2010,
+ note = {{ISBN} 3-900051-07-0},
+ url = {http://www.R-project.org/},
+}
+
+ at Manual{CRAN:inline,
+ title = {inline: Inline C, C++, Fortran function calls from
+ R},
+ author = {Oleg Sklyar and Duncan Murdoch and Mike Smith and
+ Dirk Eddelbuettel and Romain Fran\c{c}ois},
+ year = 2010,
+ note = {R package version 0.3.7},
+ url = CRAN # "/package=inline"
+}
+
+ at Manual{CRAN:Rcpp,
+ title = {Rcpp {R/C++} interface package},
+ author = {Dirk Eddelbuettel and Romain Fran\c{c}ois},
+ year = {2010},
+ note = {R package version 0.8.8},
+ url = {http://CRAN.R-project.org/package=Rcpp}
+}
+
+ at Manual{CRAN:RcppArmadillo,
+ title = {RcppArmadillo: Rcpp integration for Armadillo
+ templated linear algebra library},
+ author = {Romain Fran\c{c}ois and Dirk Eddelbuettel and
+ Douglas Bates},
+ year = 2010,
+ note = {R package version 0.2.9},
+ url = CRAN # "/package=RcppArmadillo"
+}
+
+ at Manual{CRAN:DEoptim,
+ title = {DEoptim: Global optimization by differential evolution},
+ author = {David Ardia and Katharine Mullen and Brian Peterson and Joshua Ulrich},
+ year = 2010,
+ note = {R package version 2.0-7},
+ url = CRAN # "/package=DEoptim"
+}
+
+ at Unpublished{MullenArdiaEtAl:2009:DEoptim,
+ author = {Katherine M. Mullen and David Ardia and David L. Gil
+ and Donald Windover and James Cline},
+ title = {{DEoptim}: An 'R' Package for Global Optimization by
+ Differential Evolution},
+ note = {Unpublished Manuscript},
+ year = 2009,
+ url = "http://ssrn.com/abstract=1526466"
+}
+
+ at Unpublished{ArdiaBoudtCarlEtAl:2010:DEoptim,
+ author = {David Ardia and Kris Boudt and Peter Carl and
+ Katherine M. Mullen and Brian G. Peterson},
+ title = {Differential Evolution (DEoptim) for Non-Convex
+ Portfolio Optimization},
+ note = {Unpublished Manuscript},
+ year = 2010,
+ url = "http://ssrn.com/abstract=1584905"
+}
+
+ at Book{PriceStornLampinen:2006:DE,
+ author = {Kenneth V. Price and Rainer M. Storn and Jouni
+ A. Lampinen},
+ title = {Differential Evolution -- A Practical Approach to
+ Global Optimization},
+ publisher = {Springer},
+ year = 2006,
+ address = {Berlin and Heidelberg},
+ note = {ISBN 3540209506}
+}
+
+ at TechReport{Sanderson:2010:Armadillo,
+ author = {Conrad Sanderson},
+ title = {{Armadillo}: {An} open source {C++} Algebra Library
+ for Fast Prototyping and Computationally Intensive
+ Experiments },
+ institution = {{NICTA}},
+ year = 2010,
+ url = "http://arma.sf.net"
+}
Added: pkg/RcppDE/inst/doc/RcppDE.tex
===================================================================
--- pkg/RcppDE/inst/doc/RcppDE.tex (rev 0)
+++ pkg/RcppDE/inst/doc/RcppDE.tex 2010-12-02 04:59:12 UTC (rev 2642)
@@ -0,0 +1,1581 @@
+%% use JSS class but for now with nojss option
+\documentclass[nojss,shortnames,article]{jss}
+\usepackage{rotating}
+\usepackage{float}
+\usepackage{booktabs}
+
+\author{Dirk Eddelbuettel\\Debian Project} % \And Second Author\\Plus Affiliation}
+\title{From \pkg{DEoptim} to \pkg{RcppDE}: \\
+ A case study in porting from \proglang{C} to \proglang{C++} \\
+ using \pkg{Rcpp} and \pkg{RcppArmadillo}}
+
+\Plainauthor{Dirk Eddelbuettel} % , Second Author} %% comma-separated
+\Plaintitle{DEoptim: A case study in porting to C++ and Rcpp}
+\Shorttitle{A case study in porting to C++ and Rcpp}
+
+\Abstract{
+ \noindent
+ \pkg{DEoptim} \citep{MullenArdiaEtAl:2009:DEoptim,ArdiaBoudtCarlEtAl:2010:DEoptim,CRAN:DEoptim}
+ provides differential evolution optimisation for
+ \proglang{R}. It is based on an implementation by Storn
+ \citep{PriceStornLampinen:2006:DE} and was originally implemented as an
+ interpreted \proglang{R} script. It was then rewritten in ANSI C which
+ resulted in a much improved performance.
+
+ The present paper introduces another implementation. This version is
+ written in \proglang{C++} based on the \pkg{Rcpp} package \citep{CRAN:Rcpp}
+ which provides tools for a more direct integration of \proglang{R} objects at the \proglang{C++}
+ level---and vice versa. It also uses the \pkg{RcppArmadillo} package
+ which provides an interface from \proglang{R} to the \pkg{Armadillo} linear algebra
+ package written in \proglang{C++} by Sanderson \citep{Sanderson:2010:Armadillo}.
+}
+
+\Keywords{\pkg{Rcpp}, \pkg{RcppArmadillo}, \pkg{DEoptim}, differential
+ evolution, genetic algorithm} %% at least one keyword must be supplied
+\Plainkeywords{Rcpp, RcppArmadillo, DEoptim, differential evolution, genetic algorith} %% without formatting
+
+%% publication information
+%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
+%% \Volume{13}
+%% \Issue{9}
+%% \Month{September}
+%% \Year{2004}
+%% \Submitdate{2004-09-29}
+%% \Acceptdate{2004-09-29}
+
+\Address{
+ Dirk Eddelbuettel \\
+ Debian Project \\
+ River Forest, IL, USA\\
+ %% Telephone: +43/1/31336-5053
+ %% Fax: +43/1/31336-734
+ E-mail: \email{edd at debian.org}\\
+ URL: \url{http://dirk.eddelbuettel.com}
+}
+
+%% need no \usepackage{Sweave.sty}
+
+% ------------------------------------------------------------------------
+
+\begin{document}
+
+\section{Introduction}
+
+\pkg{DEoptim}
+\citep{MullenArdiaEtAl:2009:DEoptim,ArdiaBoudtCarlEtAl:2010:DEoptim,CRAN:DEoptim}
+provides differential evolution optimisation for the \proglang{R} language
+and statistical environement. Differential optimisation is one of several
+evolutionary computing approaches; genetic algorithns and simulated annealing
+are two other ones. Differential optimisation is reasonably close to genetic
+algorithms but differs in one key aspect: parameter values are encoded as
+floating point values (rather than sequences of binary digits) which makes it
+particular suitable for real-valued optimisation problems.
+
+\pkg{DEoptim} is based on an implementation by Storn
+\citep{PriceStornLampinen:2006:DE} and was originally implemented as an
+interpreted \proglang{R} script. It was later rewritten in ANSI C which
+resulted in a much improved performance. \pkg{DEoptim} is being to optimise
+problems from a wide range of problem domains and is also being by two other
+packages.
+
+The present paper introduces the \proglang{R} package \pkg{RcppDE} which
+provides another implementation of differential evolution. This version is
+based closely on \pkg{DEoptim} but written in \proglang{C++} based on the
+\pkg{Rcpp} package \citep{CRAN:Rcpp} which provides tools for a more direct
+integration of \proglang{R} objects at the \proglang{C++} level---and vice
+versa. It also uses the \pkg{RcppArmadillo} package which provides an
+interface from \proglang{R} to the \pkg{Armadillo} linear algebra package
+written in \proglang{C++} by Sanderson \citep{Sanderson:2010:Armadillo}.
+
+The code structure descends directly from the current \pkg{DEoptim} by
+\cite{CRAN:DEoptim}. The conversion to \proglang{C++} was undertaken to see
+whether one or more of the goals \textsl{shorter}, \textsl{easier} and
+\textsl{faster} could be achieved by switching the implementation
+language. These goals were loosely defined as follows:
+\begin{itemize}
+\item[shorter] replacing code that is by necessity somewhat verbose when
+ written in \proglang{C} with more compact code written in \proglang{C++}:
+ an example would be copying of a matrix which is implemented as a dual loop
+ copying each element --- whereas \proglang{C++} allows us to use a single
+ (overloaded) \verb|+| operator and hence a single statement;
+\item[easier] this is a corollary to the previous point but also covers
+ aspect as the automatic type conversion offered by \pkg{Rcpp} as well as
+ the automatic memory management: by replacing allocation and freeing of
+ heap-based dynamic memory, a conistent source of programmer error would be
+ elimnated;
+\item[faster] this was a bit more of a conjecture as ultimately,
+ \proglang{C++} and \proglang{C} can be expected to be roughly equivalent
+ given matching compiler versions etc; however gains maybe be expected from
+ replacing a copying operation of block of adjacent memory cells with a
+ single \verb|memcpy()| call done behind the scenes; \pkg{RcppArmadillo}
+ also offers further possible gains from template metaprogramming which can
+ result in the elimination of temporary object in complex expression where,
+ loosely speaking, compile-time effort is substituted to gain later run-time
+ performance.
+\end{itemize}
+
+This paper is organised as follows. The next sections describes the structure
+of \pkg{DEoptim} which \pkg{RcppDE} shadows closesly. The following twp
+section compares differences at the \proglang{R} and \proglang{C++} level,
+respectively. Next, we changes in auxiliarry files are discusses before a
+short section reviews perfomance changes. A summary concludes. The
+appendix contains a list of figures contrasting the two implementations.
+
+\section[DEoptim structure]{\pkg{DEoptim} structure}
+
+\pkg{DEoptim} is straightforward and well-implemented package. Its
+functionality is provided by three \proglang{R} files, as well as three
+\proglang{C} files.
+
+Table~\ref{tab:Rfiles} lists the files and corresponding key functions:
+
+\begin{table}[htb]
+ \begin{center}
+ \begin{tabular}{lr}
+ \toprule
+ File & Function \\
+% \cmidrule{2}
+ \verb|DEoptim.R| \phantom{XXXXXXXXX} & \verb|DEoptim()| \\
+ & \verb|DEoptim.control()| \\[6pt]
+ \verb|methods.R| & \verb|summary.DEoptim()| \\
+ & \verb|plot.DEoptim()| \\[6pt]
+ \verb|zzz.R| & \verb|.onLoad()| \\
+ \bottomrule
+ \end{tabular}
+ \caption{Source file organisation for \proglang{R} files in \pkg{DEoptim}}
+ \label{tab:Rfiles}
+ \end{center}
+\end{table}
+
+Very few changes has to made for \pkg{RcppDE} as keeping the interface
+compatibale was an important goal. Section~\ref{sec:Rchanges} below discusses
+this in more detail.
+
+Similarly, table~\ref{tab:Cfiles} lists the \proglang{C} files.
+
+\begin{table}[htb]
+ \begin{center}
+ \begin{tabular}{lr}
+ \toprule
+ File & Function \\
+% \cmidrule{2}
+ \verb|de4_0.c| & \verb|DEoptimC()| \\
+ & \verb|devol()| \\
+ & \verb|permute()| \\[6pt]
+ \verb|evaluate.c| & \verb|evaluate()| \\[6pt]
+ \verb|get_element.c|\phantom{XXXXXXXXX} & \verb|getListElement()| \\
+
+ \bottomrule
+ \end{tabular}
+ \caption{Source file organisation for \proglang{C} files in \pkg{DEoptim}}
+ \label{tab:Cfiles}
+ \end{center}
+\end{table}
+
+A number of changes were made in the port to \proglang{C++},
+section~\ref{sec:Cppchanges} discusses these in more detail.
+
+\section[R changes]{\proglang{R} changes}
+\label{sec:Rchanges}
+
+TBD
+
+\section[C++ changes]{\proglang{C++} changes}
+\label{sec:Cppchanges}
+
+\subsection[de4_0.c and deoptim.cpp]{\code{de4\_0.c} and \code{deoptim.cpp}}
+
+\subsection[de4_0.c and devol.cpp]{\code{de4\_0.c} and \code{devol.cpp}}
+
+\bibliography{RcppDE}
+
+
+\section*{Appendix}
+
+\begin{sidewaysfigure} % fig 1: beginning of DEoptimC / DEoptim
+ \begin{minipage}{0.40\linewidth}
+ \tiny
+ \begin{CodeChunk}
+ \begin{CodeInput}
+SEXP DEoptimC(SEXP lower, SEXP upper, SEXP fn, SEXP control, SEXP rho)
+{
+ int i, j;
+
+ /* External pointers to return to R */
+ SEXP sexp_bestmem, sexp_bestval, sexp_nfeval, sexp_iter,
+ out, out_names, sexp_pop, sexp_storepop, sexp_bestmemit, sexp_bestvalit;
+
+ if (!isFunction(fn))
+ error("fn is not a function!");
+ if (!isEnvironment(rho))
+ error("rho is not an environment!");
+
+ /*-----Initialization of annealing parameters-------------------------*/
+ /* value to reach */
+ double VTR = NUMERIC_VALUE(getListElement(control, "VTR"));
+ /* chooses DE-strategy */
+ int i_strategy = INTEGER_VALUE(getListElement(control, "strategy"));
+ /* Maximum number of generations */
+ int i_itermax = INTEGER_VALUE(getListElement(control, "itermax"));
+ /* Number of objective function evaluations */
+ long l_nfeval = (long)NUMERIC_VALUE(getListElement(control, "nfeval"));
+ /* Dimension of parameter vector */
+ int i_D = INTEGER_VALUE(getListElement(control, "npar"));
+ /* Number of population members */
+ int i_NP = INTEGER_VALUE(getListElement(control, "NP"));
+ /* When to start storing populations */
+ int i_storepopfrom = INTEGER_VALUE(getListElement(control, "storepopfrom"))-1;
+ /* How often to store populations */
+ int i_storepopfreq = INTEGER_VALUE(getListElement(control, "storepopfreq"));
+ /* User-defined inital population */
+ int i_specinitialpop = INTEGER_VALUE(getListElement(control, "specinitialpop"));
+ double *initialpopv = NUMERIC_POINTER(getListElement(control, "initialpop"));
+ /* User-defined bounds */
+ double *f_lower = NUMERIC_POINTER(lower);
+ double *f_upper = NUMERIC_POINTER(upper);
+ /* stepsize */
+ double f_weight = NUMERIC_VALUE(getListElement(control, "F"));
+ /* crossover probability */
+ double f_cross = NUMERIC_VALUE(getListElement(control, "CR"));
+ /* Best of parent and child */
+ int i_bs_flag = NUMERIC_VALUE(getListElement(control, "bs"));
+ /* Print progress? */
+ int i_trace = NUMERIC_VALUE(getListElement(control, "trace"));
+ /* Re-evaluate best parameter vector? */
+ int i_check_winner = NUMERIC_VALUE(getListElement(control, "checkWinner"));
+ /* Average */
+ int i_av_winner = NUMERIC_VALUE(getListElement(control, "avWinner"));
+ /* p to define the top 100p% best solutions */
+ double i_pPct = NUMERIC_VALUE(getListElement(control, "p"));
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel A: \proglang{C} version}
+ \tiny
+
+ \end{minipage}
+ \begin{minipage}{0.03\linewidth}
+ \phantom{XX}
+ \end{minipage}
+ \begin{minipage}{0.56\linewidth}
+ \tiny
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+RcppExport SEXP DEoptim(SEXP lowerS, SEXP upperS, SEXP fnS, SEXP controlS, SEXP rhoS) {
+
+ try {
+ Rcpp::NumericVector f_lower(lowerS), f_upper(upperS); // User-defined bounds
+ Rcpp::List control(controlS); // named list of params
+
+ double VTR = Rcpp::as<double>(control["VTR"]); // value to reach
+ int i_strategy = Rcpp::as<int>(control["strategy"]); // chooses DE-strategy
+ int i_itermax = Rcpp::as<int>(control["itermax"]); // Maximum number of generations
+ long l_nfeval = 0; // nb of function evals (NOT passed in)
+ int i_D = Rcpp::as<int>(control["npar"]); // Dimension of parameter vector
+ int i_NP = Rcpp::as<int>(control["NP"]); // Number of population members
+ int i_storepopfrom = Rcpp::as<int>(control["storepopfrom"]) - 1; // When to start storing populations
+ int i_storepopfreq = Rcpp::as<int>(control["storepopfreq"]); // How often to store populations
+ int i_specinitialpop = Rcpp::as<int>(control["specinitialpop"]);// User-defined inital population
+ Rcpp::NumericMatrix initialpopm = Rcpp::as<Rcpp::NumericMatrix>(control["initialpop"]);
+ double f_weight = Rcpp::as<double>(control["F"]); // stepsize
+ double f_cross = Rcpp::as<double>(control["CR"]); // crossover probability
+ int i_bs_flag = Rcpp::as<int>(control["bs"]); // Best of parent and child
+ int i_trace = Rcpp::as<int>(control["trace"]); // Print progress?
+ int i_check_winner = Rcpp::as<int>(control["checkWinner"]); // Re-evaluate best parameter vector?
+ int i_av_winner = Rcpp::as<int>(control["avWinner"]); // Average
+ double i_pPct = Rcpp::as<double>(control["p"]); // p to define the top 100p% best solutions
+
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel B: \proglang{C++} version using \pkg{Rcpp}}
+ \end{minipage}
+ \caption{Beginning of \code{DEoptim()} \proglang{C/C++} function}
+ \label{fig:deoptim_start1}
+\end{sidewaysfigure}
+
+\begin{sidewaysfigure} % fig 2: memory allocations
+ \begin{minipage}{0.40\linewidth}
+ \tiny
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ /* Data structures for parameter vectors */
+ double **gta_popP = (double **)R_alloc(i_NP*2,sizeof(double *));
+ for (int i = 0; i < (i_NP*2); i++)
+ gta_popP[i] = (double *)R_alloc(i_D,sizeof(double));
+
+ double **gta_oldP = (double **)R_alloc(i_NP,sizeof(double *));
+ for (int i = 0; i < i_NP; i++)
+ gta_oldP[i] = (double *)R_alloc(i_D,sizeof(double));
+
+ double **gta_newP = (double **)R_alloc(i_NP,sizeof(double *));
+ for (int i = 0; i < i_NP; i++)
+ gta_newP[i] = (double *)R_alloc(i_D,sizeof(double));
+
+ double *gt_bestP = (double *)R_alloc(1,sizeof(double) * i_D);
+
+ /* Data structures for objective function values associated with
+ * parameter vectors */
+ double *gta_popC = (double *)R_alloc(i_NP*2,sizeof(double));
+ double *gta_oldC = (double *)R_alloc(i_NP,sizeof(double));
+ double *gta_newC = (double *)R_alloc(i_NP,sizeof(double));
+ double *gt_bestC = (double *)R_alloc(1,sizeof(double));
+
+ double *t_bestitP = (double *)R_alloc(1,sizeof(double) * i_D);
+ double *t_tmpP = (double *)R_alloc(1,sizeof(double) * i_D);
+ double *tempP = (double *)R_alloc(1,sizeof(double) * i_D);
+
+ int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
+ double *gd_pop = (double *)R_alloc(i_NP*i_D,sizeof(double));
+ double *gd_storepop = (double *)R_alloc(i_NP,sizeof(double) * i_D * i_nstorepop);
+ double *gd_bestmemit = (double *)R_alloc(i_itermax*i_D,sizeof(double));
+ double *gd_bestvalit = (double *)R_alloc(i_itermax,sizeof(double));
+ int gi_iter = 0;
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel A: \proglang{C} version}
+ \tiny
+ \end{minipage}
+ \begin{minipage}{0.03\linewidth}
+ \phantom{XX}
+ \end{minipage}
+ \begin{minipage}{0.56\linewidth}
+ \tiny
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ arma::colvec minbound(f_lower.begin(), f_lower.size(), false); // convert Rcpp vectors to arma vectors
+ arma::colvec maxbound(f_upper.begin(), f_upper.size(), false);
+ arma::mat initpopm(initialpopm.begin(), initialpopm.rows(), initialpopm.cols(), false);
+
+ arma::mat ta_popP(i_D, i_NP*2); // Data structures for parameter vectors
+ arma::mat ta_oldP(i_D, i_NP);
+ arma::mat ta_newP(i_D, i_NP);
+ arma::colvec t_bestP(i_D);
+
+ arma::colvec ta_popC(i_NP*2); // Data structures for obj. fun. values
+ arma::colvec ta_oldC(i_NP);
+ arma::colvec ta_newC(i_NP);
+ double t_bestC;
+
+ arma::colvec t_bestitP(i_D);
+ arma::colvec t_tmpP(i_D);
+
+ int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
+ arma::mat d_pop(i_D, i_NP);
+ Rcpp::List d_storepop(i_nstorepop);
+ arma::mat d_bestmemit(i_D, i_itermax);
+ arma::colvec d_bestvalit(i_itermax);
+ int i_iter = 0;
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel B: \proglang{C++} version using \pkg{Rcpp}}
+ \end{minipage}
+ \caption{Memory allocation in \code{DEoptim()} \proglang{C/C++} function}
+ \label{fig:deoptim_memory}
+\end{sidewaysfigure}
+
+
+\begin{sidewaysfigure} % fig 3: devol and arguments
+ \begin{minipage}{0.40\linewidth}
+ \tiny
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ /*---optimization--------------------------------------*/
+ devol(VTR, f_weight, f_cross, i_bs_flag, f_lower, f_upper, fn, rho, i_trace,
+ i_strategy, i_D, i_NP, i_itermax,
+ initialpopv, i_storepopfrom, i_storepopfreq,
+ i_specinitialpop, i_check_winner, i_av_winner,
+ gta_popP, gta_oldP, gta_newP, gt_bestP,
+ gta_popC, gta_oldC, gta_newC, gt_bestC,
+ t_bestitP, t_tmpP, tempP,
+ gd_pop, gd_storepop, gd_bestmemit, gd_bestvalit,
+ &gi_iter, i_pPct, &l_nfeval);
+ /*---end optimization----------------------------------*/
+
+ PROTECT(sexp_bestmem = NEW_NUMERIC(i_D));
+ for (i = 0; i < i_D; i++) {
+ NUMERIC_POINTER(sexp_bestmem)[i] = gt_bestP[i];
+ }
+
+ j = i_NP * i_D;
+ PROTECT(sexp_pop = NEW_NUMERIC(j));
+ for (i = 0; i < j; i++)
+ NUMERIC_POINTER(sexp_pop)[i] = gd_pop[i];
+
+ j = i_nstorepop * i_NP * i_D;
+ PROTECT(sexp_storepop = NEW_NUMERIC(j));
+ for (i = 0; i < j; i++)
+ NUMERIC_POINTER(sexp_storepop)[i] = gd_storepop[i];
+
+ j = gi_iter * i_D;
+ PROTECT(sexp_bestmemit = NEW_NUMERIC(j));
+ for (i = 0; i < j; i++)
+ NUMERIC_POINTER(sexp_bestmemit)[i] = gd_bestmemit[i];
+ j = gi_iter;
+ PROTECT(sexp_bestvalit = NEW_NUMERIC(j));
+ for (i = 0; i < j; i++)
+ NUMERIC_POINTER(sexp_bestvalit)[i] = gd_bestvalit[i];
+
+ PROTECT(sexp_bestval = NEW_NUMERIC(1));
+ NUMERIC_POINTER(sexp_bestval)[0] = gt_bestC[0];
+
+ PROTECT(sexp_nfeval = NEW_INTEGER(1));
+ //INTEGER_POINTER(sexp_nfeval)[0] = 0;
+ INTEGER_POINTER(sexp_nfeval)[0] = l_nfeval;
+
+ PROTECT(sexp_iter = NEW_INTEGER(1));
+ INTEGER_POINTER(sexp_iter)[0] = gi_iter;
+
+ PROTECT(out = NEW_LIST(8));
+ SET_VECTOR_ELT(out, 0, sexp_bestmem);
+ SET_VECTOR_ELT(out, 1, sexp_bestval);
+ SET_VECTOR_ELT(out, 2, sexp_nfeval);
+ SET_VECTOR_ELT(out, 3, sexp_iter);
+ SET_VECTOR_ELT(out, 4, sexp_bestmemit);
+ SET_VECTOR_ELT(out, 5, sexp_bestvalit);
+ SET_VECTOR_ELT(out, 6, sexp_pop);
+ SET_VECTOR_ELT(out, 7, sexp_storepop);
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel A: \proglang{C} version}
+ \tiny
+ \end{minipage}
+ \begin{minipage}{0.03\linewidth}
+ \phantom{XX}
+ \end{minipage}
+ \begin{minipage}{0.56\linewidth}
+ \tiny
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ PROTECT(out_names = NEW_STRING(8));
+ SET_STRING_ELT(out_names, 0, mkChar("bestmem"));
+ SET_STRING_ELT(out_names, 1, mkChar("bestval"));
+ SET_STRING_ELT(out_names, 2, mkChar("nfeval"));
+ SET_STRING_ELT(out_names, 3, mkChar("iter"));
+ SET_STRING_ELT(out_names, 4, mkChar("bestmemit"));
+ SET_STRING_ELT(out_names, 5, mkChar("bestvalit"));
+ SET_STRING_ELT(out_names, 6, mkChar("pop"));
+ SET_STRING_ELT(out_names, 7, mkChar("storepop"));
+
+ SET_NAMES(out, out_names);
+
+ UNPROTECT(10);
+
+ return out;
+}
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel A: \proglang{C} version (in both columns)}
+ \tiny
+
+ \bigskip
+ \phantom{XX}
+ \bigskip
+ \phantom{XX}
+ \bigskip
+ \phantom{XX}
+ \bigskip
+
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ // call actual Differential Evolution optimization given the parameters
+ devol(VTR, f_weight, f_cross, i_bs_flag, minbound, maxbound, fnS, rhoS, i_trace, i_strategy, i_D, i_NP,
+ i_itermax, initpopm, i_storepopfrom, i_storepopfreq, i_specinitialpop, i_check_winner, i_av_winner,
+ ta_popP, ta_oldP, ta_newP, t_bestP, ta_popC, ta_oldC, ta_newC, t_bestC, t_bestitP, t_tmpP,
+ d_pop, d_storepop, d_bestmemit, d_bestvalit, i_iter, i_pPct, l_nfeval);
+
+ return Rcpp::List::create(Rcpp::Named("bestmem") = t_bestP, // and return a named list with results to R
+ Rcpp::Named("bestval") = t_bestC,
+ Rcpp::Named("nfeval") = l_nfeval,
+ Rcpp::Named("iter") = i_iter,
+ Rcpp::Named("bestmemit") = trans(d_bestmemit),
+ Rcpp::Named("bestvalit") = d_bestvalit,
+ Rcpp::Named("pop") = trans(d_pop),
+ Rcpp::Named("storepop") = d_storepop);
+
+ } catch( std::exception& ex) {
+ forward_exception_to_r(ex);
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)");
+ }
+ return R_NilValue;
+}
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel B: \proglang{C++} version using \pkg{Rcpp}}
+ \end{minipage}
+ \caption{\code{DEoptim()} call of \code{devol()} and return of results to \proglang{R}}
+ \label{fig:deoptim_end}
+\end{sidewaysfigure}
+
+
+
+\begin{sidewaysfigure} % fig 4: devol and arguments
+ \begin{minipage}{0.40\linewidth}
+ \tiny
+ \begin{CodeChunk}
+ \begin{CodeInput}
+void devol(double VTR, double f_weight, double f_cross, int i_bs_flag,
+ double *lower, double *upper, SEXP fcall, SEXP rho, int trace,
+ int i_strategy, int i_D, int i_NP, int i_itermax,
+ double *initialpopv, int i_storepopfrom, int i_storepopfreq,
+ int i_specinitialpop, int i_check_winner, int i_av_winner,
+ double **gta_popP, double **gta_oldP, double **gta_newP, double *gt_bestP,
+ double *gta_popC, double *gta_oldC, double *gta_newC, double *gt_bestC,
+ double *t_bestitP, double *t_tmpP, double *tempP,
+ double *gd_pop, double *gd_storepop, double *gd_bestmemit, double *gd_bestvalit,
+ int *gi_iter, double i_pPct, long *l_nfeval)
+{
+
+#define URN_DEPTH 5 /* 4 + one index to avoid */
+
+ /* initialize parameter vector to pass to evaluate function */
+ SEXP par; PROTECT(par = NEW_NUMERIC(i_D));
+
+ int i, j, k, x; /* counting variables */
+ int i_r1, i_r2, i_r3, i_r4; /* placeholders for random indexes */
+
+ int ia_urn2[URN_DEPTH];
+ int i_nstorepop, i_xav;
+ i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
+
+ int popcnt, bestacnt, same; /* lazy cnters */
+
+ double *fa_minbound = lower;
+ double *fa_maxbound = upper;
+ double f_jitter, f_dither;
+
+ double t_bestitC;
+ double t_tmpC, tmp_best;
+
+ double initialpop[i_NP][i_D];
+
+ /* vars for DE/current-to-p-best/1 */
+ int i_pbest;
+ int p_NP = round(i_pPct * i_NP); /* choose at least two best solutions */
+ p_NP = p_NP < 2 ? 2 : p_NP;
+ int sortIndex[i_NP]; /* sorted values of gta_oldC */
+ for(i = 0; i < i_NP; i++) sortIndex[i] = i;
+
+ /* vars for when i_bs_flag == 1 */
+ int i_len, done, step, bound;
+ double tempC;
+
+ GetRNGstate();
+
+ gta_popP[0][0] = 0;
+ \end{CodeInput}
+ \end{CodeChunk}
+ \normalsize
+ \centering{Panel A: \proglang{C} version}
+ \tiny
+ \end{minipage}
+ \begin{minipage}{0.03\linewidth}
+ \phantom{XX}
+ \end{minipage}
+ \begin{minipage}{0.56\linewidth}
+ \tiny
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+void devol(double VTR, double f_weight, double f_cross, int i_bs_flag,
+ arma::colvec & fa_minbound, arma::colvec & fa_maxbound, SEXP fcall, SEXP rho, int i_trace,
+ int i_strategy, int i_D, int i_NP, int i_itermax, arma::mat & initialpopm,
+ int i_storepopfrom, int i_storepopfreq, int i_specinitialpop, int i_check_winner, int i_av_winner,
+ arma::mat &ta_popP, arma::mat &ta_oldP, arma::mat &ta_newP, arma::colvec & t_bestP,
+ arma::colvec & ta_popC, arma::colvec & ta_oldC, arma::colvec & ta_newC, double & t_bestC,
+ arma::colvec & t_bestitP, arma::colvec & t_tmpP,
+ arma::mat &d_pop, Rcpp::List &d_storepop, arma::mat & d_bestmemit, arma::colvec & d_bestvalit,
+ int & i_iterations, double i_pPct, long & l_nfeval) {
+
+ Rcpp::DE::EvalBase *ev = NULL; // pointer to abstract base class
+ if (TYPEOF(fcall) == EXTPTRSXP) { // non-standard mode: we are being passed an external pointer
+ ev = new Rcpp::DE::EvalCompiled(fcall); // so assign a pointer using external pointer in fcall SEXP
+ } else { // standard mode: env_ is an env, fcall_ is a function
+ ev = new Rcpp::DE::EvalStandard(fcall, rho); // so assign R function and environment
+ }
+ const int urn_depth = 5; // 4 + one index to avoid
+ Rcpp::NumericVector par(i_D); // initialize parameter vector to pass to evaluate function
+ arma::icolvec::fixed<urn_depth> ia_urn2; // fixed-size vector for urn draws
+ arma::icolvec ia_urntmp(i_NP); // so that we don't need to re-allocated each time in permute
+ arma::mat initialpop(i_D, i_NP);
+ int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
+ int p_NP = round(i_pPct * i_NP); // choose at least two best solutions
+ p_NP = p_NP < 2 ? 2 : p_NP;
+ arma::icolvec sortIndex(i_NP); // sorted values of ta_oldC
+ if (i_strategy == 6) {
+ for (int i = 0; i < i_NP; i++)
+ sortIndex[i] = i;
+ }
+ GetRNGstate();
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel B: \proglang{C++} version using \pkg{Rcpp}}
+ \end{minipage}
+ \caption{\code{devol()} beginning}
+ \label{fig:deoptim_end}
+\end{sidewaysfigure}
+
+
+
+\begin{sidewaysfigure} % fig 5: devol inits
+ \begin{minipage}{0.40\linewidth}
+ \tiny
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ /* initialize initial popuplation */
+ for (int i = 0; i < i_NP; i++) {
+ for (int j = 0; j < i_D; j++) {
+ initialpop[i][j] = 0.0;
+ }
+ }
+
+ /* initialize best members */
+ for (int i = 0; i < i_itermax * i_D; i++)
+ gd_bestmemit[i] = 0.0;
+
+ /* initialize best values */
+ for (int i = 0; i < i_itermax; i++)
+ gd_bestvalit[i] = 0.0;
+
+ /* initialize best population */
+ for (int i = 0; i < i_NP * i_D; i++)
+ gd_pop[i] = 0.0;
+
+ /* initialize stored populations */
+ if (i_nstorepop < 0)
+ i_nstorepop = 0;
+
+ for (int i = 0; i < (i_nstorepop * i_NP * i_D); i++)
+ gd_storepop[i] = 0.0;
+
+ /* if initial population provided, initialize with values */
+ if (i_specinitialpop > 0) {
+ k = 0;
+
+ for (j = 0; j < i_D; j++) {
+ for (i = 0; i < i_NP; i++) {
+ initialpop[i][j] = initialpopv[k];
+ k += 1;
+ }
+ }
+ }
+
+ /* number of function evaluations
+ * (this is an input via DEoptim.control, but we over-write it?) */
+ *l_nfeval = 0;
+
+ /*------Initialization-----------------------------*/
+ for (i = 0; i < i_NP; i++) {
+ for (j = 0; j < i_D; j++) {
+ if (i_specinitialpop <= 0) { /* random initial member */
+ gta_popP[i][j] = fa_minbound[j] +
+ unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
+
+ }
+ else /* or user-specified initial member */
+ gta_popP[i][j] = initialpop[i][j];
+ }
+ gta_popC[i] = evaluate(l_nfeval, gta_popP[i], par, fcall, rho);
+
+ if (i == 0 || gta_popC[i] <= gt_bestC[0]) {
+ gt_bestC[0] = gta_popC[i];
+ for (j = 0; j < i_D; j++)
+ gt_bestP[j]=gta_popP[i][j];
+ }
+ }
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ %\normalsize \centering{Panel A: \proglang{C} version} \tiny
+ \end{minipage}
+ \begin{minipage}{0.03\linewidth}
+ \phantom{XX}
+ \end{minipage}
+ \begin{minipage}{0.56\linewidth}
+ \tiny
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+
+ /*---assign pointers to current ("old") population---*/
+ gta_oldP = gta_popP;
+ gta_oldC = gta_popC;
+
+ /*------Iteration loop--------------------------------------------*/
+ int i_iter = 0;
+ popcnt = 0;
+ bestacnt = 0;
+ i_xav = 1;
+ \end{CodeInput}
+ \end{CodeChunk}
+
+ \normalsize
+ \centering{Panel A: \proglang{C} version (in both columns)}
+ \tiny
+
+ \bigskip
+ \phantom{XX}
+ \bigskip
+ \phantom{XX}
+ \bigskip
+ \phantom{XX}
+ \bigskip
+
+
+ \begin{CodeChunk}
+ \begin{CodeInput}
+ initialpop.zeros(); // initialize initial popuplation
+ d_bestmemit.zeros(); // initialize best members
+ d_bestvalit.zeros(); // initialize best values
+ d_pop.zeros(); // initialize best population
+ i_nstorepop = (i_nstorepop < 0) ? 0 : i_nstorepop;
+
+ if (i_specinitialpop > 0) { // if initial population provided, initialize with values
+ initialpop = trans(initialpopm); // transpose as we prefer columns for population members here
+ }
+
+ for (int i = 0; i < i_NP; i++) { // ------Initialization-----------------------------
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 2642
More information about the Rcpp-commits
mailing list