[Rcpp-commits] r2430 - in pkg: . RcppDE/scripts RcppDE/src RcppParDE RcppParDE/R RcppParDE/man RcppParDE/scripts RcppParDE/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 12 04:28:42 CET 2010
Author: edd
Date: 2010-11-12 04:28:41 +0100 (Fri, 12 Nov 2010)
New Revision: 2430
Added:
pkg/RcppParDE/
pkg/RcppParDE/DESCRIPTION
pkg/RcppParDE/NAMESPACE
pkg/RcppParDE/R/
pkg/RcppParDE/R/DEoptim.R
pkg/RcppParDE/R/methods.R
pkg/RcppParDE/R/zzz.R
pkg/RcppParDE/man/
pkg/RcppParDE/man/DEoptim-methods.Rd
pkg/RcppParDE/man/DEoptim.Rd
pkg/RcppParDE/man/DEoptim.control.Rd
pkg/RcppParDE/scripts/
pkg/RcppParDE/scripts/openmp.r
pkg/RcppParDE/src/
pkg/RcppParDE/src/Makevars
pkg/RcppParDE/src/deoptim.cpp
pkg/RcppParDE/src/devol.cpp
pkg/RcppParDE/src/evaluate.cpp
pkg/RcppParDE/src/evaluate.h
pkg/RcppParDE/src/permute.cpp
Removed:
pkg/RcppDE/scripts/openmp.r
pkg/RcppDE/src/devolMP.cpp
pkg/RcppDE/src/evaluateMP.cpp
pkg/RcppDE/src/permuteMP.cpp
Modified:
pkg/RcppDE/src/Makevars
pkg/RcppDE/src/deoptim.cpp
pkg/RcppDE/src/devol.cpp
pkg/RcppDE/src/permute.cpp
Log:
splitting RcppParDE off RcppDE so that we can compare the approach with OpenMP (in RcppParDE) to the plain C++ approach (in RcppDE)
Deleted: pkg/RcppDE/scripts/openmp.r
===================================================================
--- pkg/RcppDE/scripts/openmp.r 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/scripts/openmp.r 2010-11-12 03:28:41 UTC (rev 2430)
@@ -1,79 +0,0 @@
-#!/usr/bin/r -t
-#
-# with OpenMP we do not get the same uniform random number draws as we act in parallel, so just compare results (and timings)
-
-suppressMessages(library(DEoptim)) # the original, currently 2.0.7
-suppressMessages(library(RcppDE)) # the contender
-suppressMessages(library(inline)) # helps with
-
-Genrose <- function(x) { ## One generalization of the Rosenbrock banana valley function (n parameters)
- n <- length(x)
- 1.0 + sum (100 * (x[-n]^2 - x[-1])^2 + (x[-1] - 1)^2)
-}
-
-
-#maxIt <- 500 # not excessive but so that we get some run-time on simple problems
-#n <- 40
-
-maxIt <- 500
-n <- 20
-
-inc <- 'double genrose(SEXP xs) {
- Rcpp::NumericVector x(xs);
- int n = x.size();
- double sum = 1.0;
- for (int i=1; i<n; i++) {
- sum += 100*( pow(x[i-1]*x[i-1] - x[i], 2)) + (x[i] - 1)*(x[i] - 1);
- }
- return(sum);
- }
-
- double wild(SEXP xs) {
- Rcpp::NumericVector x(xs);
- int n = x.size();
- double sum = 0.0;
- for (int i=0; i<n; i++) {
- sum += 10 * sin(0.3 * x[i]) * sin(1.3 * x[i]*x[i]) + 0.00001 * x[i]*x[i]*x[i]*x[i] + 0.2 * x[i] + 80;
- }
- sum /= n;
- return(sum);
- }
-
- double rastrigin(SEXP xs) {
- Rcpp::NumericVector x(xs);
- int n = x.size();
- double sum = 20.0;
- for (int i=0; i<n; i++) {
- sum += x[i]+2 - 10*cos(2*M_PI*x[i]);
- }
- return(sum);
- }
-
- '
-
-## now via a class returning external pointer
-src.xptr <- 'std::string fstr = Rcpp::as<std::string>(funname);
- typedef double (*funcPtr)(SEXP);
- if (fstr == "genrose")
- return(XPtr<funcPtr>(new funcPtr(&genrose)));
- else if (fstr == "wild")
- return(XPtr<funcPtr>(new funcPtr(&wild)));
- else
- return(XPtr<funcPtr>(new funcPtr(&rastrigin)));
- '
-
-create_xptr <- cxxfunction(signature(funname="character"), body=src.xptr, inc=inc, plugin="Rcpp")
-xptr <- create_xptr("genrose")
-
-set.seed(42)
-print(system.time( {
- res <- RcppDE::DEoptim(fn=xptr, lower=rep(-25, n), upper=rep(25, n), control=list(NP=10*n, itermax=maxIt, trace=FALSE))
- print(res[[1]][["bestval"]])
-}))
-
-set.seed(42)
-print(system.time( {
- res <- DEoptim::DEoptim(fn=Genrose, lower=rep(-25, n), upper=rep(25, n), control=list(NP=10*n, itermax=maxIt, trace=FALSE))
- print(res[[1]][["bestval"]])
-}))
-
Modified: pkg/RcppDE/src/Makevars
===================================================================
--- pkg/RcppDE/src/Makevars 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/Makevars 2010-11-12 03:28:41 UTC (rev 2430)
@@ -1,11 +1,7 @@
## Hey Emacs make this a -*- mode: makefile; -*- file
##
-## -- for OpenMP (with -D macro to switch to OpenMP enabled source file)
-PKG_CXXFLAGS=-fopenmp -DUSE_OPENMP
-PKG_LIBS= -fopenmp -lgomp $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
-##
## -- for Google Perftools profiling
## PKG_LIBS= $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lprofiler
##
## -- default
-#PKG_LIBS= $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS= $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
Modified: pkg/RcppDE/src/deoptim.cpp
===================================================================
--- pkg/RcppDE/src/deoptim.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/deoptim.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -29,7 +29,7 @@
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 evaluations (NOT passed in)
+ 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
@@ -44,7 +44,7 @@
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
- arma::colvec minbound(f_lower.begin(), f_lower.size(), false); // convert three Rcpp vectors to arma vectors
+ 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);
Modified: pkg/RcppDE/src/devol.cpp
===================================================================
--- pkg/RcppDE/src/devol.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/devol.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -7,10 +7,9 @@
// and based on DE-Engine v4.0, Rainer Storn, 2004
// (http://www.icsi.berkeley.edu/~storn/DeWin.zip)
-#ifndef USE_OPENMP
#include <RcppArmadillo.h> // declarations for both Rcpp and RcppArmadillo offering Armadillo classes
-#include <omp.h> // OpenMP for compiler-generated multithreading
#include "evaluate.h" // simple function evaluation framework
+// #include <google/profiler.h>
void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urntmp[]);
@@ -273,4 +272,3 @@
PutRNGstate();
// ProfilerStop();
}
-#endif
Deleted: pkg/RcppDE/src/devolMP.cpp
===================================================================
--- pkg/RcppDE/src/devolMP.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/devolMP.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -1,298 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
-//
-// Port of DEoptim (2.0.7) by Ardia et al to Rcpp/RcppArmadillo/Armadillo
-// Copyright (C) 2010 Dirk Eddelbuettel <edd at debian.org>
-//
-// DEoptim is Copyright (C) 2009 David Ardia and Katharine Mullen
-// and based on DE-Engine v4.0, Rainer Storn, 2004
-// (http://www.icsi.berkeley.edu/~storn/DeWin.zip)
-
-#ifdef USE_OPENMP
-#include <RcppArmadillo.h> // declarations for both Rcpp and RcppArmadillo offering Armadillo classes
-#include <omp.h> // OpenMP for compiler-generated multithreading
-#include "evaluate.h" // simple function evaluation framework
-
-void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urntmp[]);
-
-inline double drndu(void) {
-#ifdef USE_OPENMP
- return arma::randu();
-#else
- return ::unif_rand();
-#endif
-}
-
-inline int irndu(const double val) {
-#ifdef USE_OPENMP
- return static_cast<int>(arma::randu() * val);
-#else
- return static_cast<int>(::unif_rand() * val);
-#endif
-}
-
-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) {
-
- //ProfilerStart("/tmp/RcppDE.prof");
- 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();
-
- 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-----------------------------
- if (i_specinitialpop <= 0) { // random initial member
- for (int j = 0; j < i_D; j++) {
- ta_popP.at(j,i) = fa_minbound[j] + drndu() * (fa_maxbound[j] - fa_minbound[j]);
- }
- } else { // or user-specified initial member
- ta_popP.col(i) = initialpop.col(i);
- }
- memcpy(REAL(par), ta_popP.colptr(i), i_D * sizeof(double));
- ta_popC[i] = ev->eval(par);
- if (i == 0 || ta_popC[i] <= t_bestC) {
- t_bestC = ta_popC[i];
- t_bestP = ta_popP.unsafe_col(i);
- }
- }
-
- ta_oldP = ta_popP.cols(0, i_NP-1); // ---assign pointers to current ("old") population---
- ta_oldC = ta_popC.rows(0, i_NP-1);
-
- int i_iter = 0; // ------Iteration loop--------------------------------------------
- int popcnt = 0;
- int i_xav = 1;
-
- for (i_iter=0; (i_iter < i_itermax) && (t_bestC > VTR); i_iter++) { // main loop =============================
- if (i_iter % i_storepopfreq == 0 && i_iter >= i_storepopfrom) { // store intermediate populations
- d_storepop[popcnt++] = Rcpp::wrap( trans(ta_oldP) );
- } // end store pop
-
- d_bestmemit.col(i_iter) = t_bestP; // store the best member
- d_bestvalit[i_iter] = t_bestC; // store the best value
- t_bestitP = t_bestP;
- //i_iter++; // increase iteration counter
-
- double f_dither = f_weight + drndu() * (1.0 - f_weight); // ----computer dithering factor ----------
-
- if (i_strategy == 6) { // ---DE/current-to-p-best/1 --------------------------------------
- arma::colvec temp_oldC = ta_oldC; // copy of ta_oldC to not change it
- rsort_with_index( temp_oldC.memptr(), sortIndex.begin(), i_NP ); // to use sortIndex later
- }
-
-#pragma omp parallel for shared(ia_urn2,ta_oldP,ta_newP,ta_newC,t_tmpP,ia_urntmp) schedule(static,1)
- for (int i = 0; i < i_NP; i++) { // ----start of loop through ensemble------------------------
-
- t_tmpP = ta_oldP.col(i); // t_tmpP is the vector to mutate and eventually select
-
- permute(ia_urn2.memptr(), urn_depth, i_NP, i, ia_urntmp.memptr()); // Pick 4 random and distinct
- int k = 0; // loop counter used in all strategies below
-
- // ===Choice of strategy=======================================================
- switch (i_strategy) {
-
- case 1: { // ---classical strategy DE/rand/1/bin-----------------------------
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(j,ia_urn2[1]) + f_weight *
- (ta_oldP.at(j,ia_urn2[2]) - ta_oldP.at(j,ia_urn2[3]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- case 2: { // ---DE/local-to-best/1/bin---------------------------------------
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target
- t_tmpP[j] = t_tmpP[j] + f_weight * (t_bestitP[j] - t_tmpP[j]) + f_weight *
- (ta_oldP.at(j,ia_urn2[2]) - ta_oldP.at(j,ia_urn2[3]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- case 3: { // ---DE/best/1/bin with jitter------------------------------------
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target
- double f_jitter = 0.0001 * drndu() + f_weight;
- t_tmpP[j] = t_bestitP[j] + f_jitter * (ta_oldP.at(j,ia_urn2[1]) - ta_oldP.at(j,ia_urn2[2]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- case 4: { // ---DE/rand/1/bin with per-vector-dither-------------------------
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target *
- t_tmpP[j] = ta_oldP.at(j,ia_urn2[1]) + (f_weight + drndu()*(1.0 - f_weight))
- * (ta_oldP.at(j,ia_urn2[2]) - ta_oldP.at(j,ia_urn2[3]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- case 5: { // ---DE/rand/1/bin with per-generation-dither---------------------
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(j,ia_urn2[1]) + f_dither
- * (ta_oldP.at(j,ia_urn2[2]) - ta_oldP.at(j,ia_urn2[3]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- case 6: { // ---DE/current-to-p-best/1 (JADE)--------------------------------
- int i_pbest = sortIndex[irndu(p_NP)]; // select from [0, 1, 2, ..., (pNP-1)]
- int j = irndu(i_D); // random parameter
- do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(j,i) + f_weight * (ta_oldP.at(j,i_pbest) - ta_oldP.at(j,i)) +
- f_weight * (ta_oldP.at(j,ia_urn2[1]) - ta_oldP.at(j,ia_urn2[2]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- break;
- }
- default: { // ---variation to DE/rand/1/bin: either-or-algorithm--------------
- int j = irndu(i_D); // random parameter
- if (drndu() < 0.5) { // differential mutation, Pmu = 0.5
- do { // add fluctuation to random target */
- t_tmpP[j] = ta_oldP.at(j,ia_urn2[1]) +
- f_weight * (ta_oldP.at(j,ia_urn2[2]) - ta_oldP.at(j,ia_urn2[3]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
-
- } else { // recombination with K = 0.5*(F+1) -. F-K-Rule
- do { // add fluctuation to random target */
- t_tmpP[j] = ta_oldP.at(j,ia_urn2[1]) + 0.5 * (f_weight + 1.0) *
- (ta_oldP.at(j,ia_urn2[2]) + ta_oldP.at(j,ia_urn2[3]) - 2 * ta_oldP.at(j,ia_urn2[1]));
- j = (j + 1) % i_D;
- } while ((drndu() < f_cross) && (++k < i_D));
- }
- break;
- }
- } // end switch (i_strategy) ...
-
- for (int j = 0; j < i_D; j++) { // boundary constr, bounce-back method was not enforcing bounds
- if (t_tmpP[j] < fa_minbound[j]) {
- t_tmpP[j] = fa_minbound[j] + drndu() * (fa_maxbound[j] - fa_minbound[j]);
- }
- if (t_tmpP[j] > fa_maxbound[j]) {
- t_tmpP[j] = fa_maxbound[j] - drndu() * (fa_maxbound[j] - fa_minbound[j]);
- }
- }
-
- // ------Trial mutation now in t_tmpP-----------------
- memcpy(REAL(par), t_tmpP.memptr(), i_D * sizeof(double));
- double t_tmpC = ev->eval(par); // Evaluate mutant in t_tmpP
- if (t_tmpC <= ta_oldC[i] || i_bs_flag) { // i_bs_flag means will choose best NP later
- ta_newP.col(i) = t_tmpP; // replace target with mutant
- ta_newC[i] = t_tmpC;
- if (t_tmpC <= t_bestC) {
- t_bestP = t_tmpP;
- t_bestC = t_tmpC;
- }
- } else {
- ta_newP.col(i) = ta_oldP.col(i);
- ta_newC[i] = ta_oldC[i];
- }
- //#ifdef USE_OPENMP
- //std::cerr << "Thread " << omp_get_thread_num() << "\n";
- //#endif
- } // End mutation loop through pop., ie the "for (i = 0; i < i_NP; i++)"
-#pragma omp barrier
-
- if (i_bs_flag) { // examine old and new pop. and take the best NP members into next generation
-
- ta_popP.cols(0, i_NP-1) = ta_oldP;
- ta_popC.rows(0, i_NP-1) = ta_oldC;
-
- ta_popP.cols(i_NP, 2*i_NP-1) = ta_newP;
- ta_popC.rows(i_NP, 2*i_NP-1) = ta_newC;
-
- int i_len = 2 * i_NP;
- int step = i_len, done; // array length
- while (step > 1) {
- step /= 2; // halve the step size
- do {
- done = 1;
- int bound = i_len - step;
- for (int j = 0; j < bound; j++) {
- int i = j + step + 1;
- if (ta_popC[j] > ta_popC[i-1]) {
- ta_popP.swap_cols(j, i-1);
- ta_popC.swap_rows(j, i-1);
- done = 0;
- } // if
- } // for
- } while (!done); // while
- } // while (step > 1)
- ta_newP = ta_popP.cols(0, i_NP-1); // now the best NP are in first NP places in gta_pop, use them
- ta_newC = ta_popC.rows(0, i_NP-1);
- } // i_bs_flag
-
- ta_oldP = ta_newP; // have selected NP mutants move on to next generation
- ta_oldC = ta_newC;
-
- if (i_check_winner) { // check if the best stayed the same, if necessary
- int same = 1;
- for (int j = 0; j < i_D; j++) {
- if (t_bestitP[j] != t_bestP[j]) {
- same = 0;
- }
- }
- if (same && i_iter > 1) {
- i_xav++;
- memcpy(REAL(par), t_bestP.memptr(), i_D * sizeof(double));
- double tmp_best = ev->eval(par);// if re-evaluation of winner
- if (i_av_winner) // possibly letting winner be the average of all past generations
- t_bestC = ((1/(double)i_xav) * t_bestC) + ((1/(double)i_xav) * tmp_best) +
- (d_bestvalit[i_iter-1] * ((double)(i_xav - 2))/(double)i_xav);
- else
- t_bestC = tmp_best;
- } else {
- i_xav = 1;
- }
- }
- t_bestitP = t_bestP;
-
- if ( (i_trace > 0) && ((i_iter % i_trace) == 0) ) {
- Rprintf("Iteration: %d bestvalit: %f bestmemit:", i_iter, t_bestC);
- for (int j = 0; j < i_D; j++)
- Rprintf("%12.6f", t_bestP[j]);
- Rprintf("\n");
- }
- } // end loop through generations
-
- d_pop = ta_oldP;
- i_iterations = i_iter;
- l_nfeval = ev->getNbEvals();
- //PutRNGstate();
- // ProfilerStop();
-}
-#endif
Deleted: pkg/RcppDE/src/evaluateMP.cpp
===================================================================
--- pkg/RcppDE/src/evaluateMP.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/evaluateMP.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -1,25 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
-//
-// Port of DEoptim (2.0.7) by Ardia et al to Rcpp/RcppArmadillo/Armadillo
-// Copyright (C) 2010 Dirk Eddelbuettel <edd at debian.org>
-//
-// DEoptim is Copyright (C) 2009 David Ardia and Katharine Mullen
-#if 0
-#ifdef USE_OPENMP
-#include <RcppArmadillo.h>
-
-// Slighly accelerated version of evaluate to evaluate function fcall with parameters param in environment env
-// Uses externally allocated par() vector into which param are copied
-double evaluate(long &l_nfeval, const double *param, SEXP par, SEXP fcall, SEXP env) {
- memcpy(REAL(par), param, Rf_nrows(par) * sizeof(double)); // -- faster: direct access _assuming_ numeric vector
- Rcpp::NumericVector x(par);
- int n = x.size();
- double sum = 1.0;
- for (int i=1; i<n; i++) {
- sum += 100*( pow(x[i-1]*x[i-1] - x[i], 2)) + (x[i] - 1)*(x[i] - 1);
- }
- l_nfeval++; // increment function evaluation count
- return(sum);
-}
-#endif
-#endif
Modified: pkg/RcppDE/src/permute.cpp
===================================================================
--- pkg/RcppDE/src/permute.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/permute.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -7,7 +7,6 @@
// and based on DE-Engine v4.0, Rainer Storn, 2004
// (http://www.icsi.berkeley.edu/~storn/DeWin.zip)
-#ifndef USE_OPENMP
#include <RcppArmadillo.h>
// Function : void permute(int ia_urn2[], int i_urn2_depth)
@@ -46,4 +45,3 @@
}
PutRNGstate();
}
-#endif
Deleted: pkg/RcppDE/src/permuteMP.cpp
===================================================================
--- pkg/RcppDE/src/permuteMP.cpp 2010-11-12 03:00:53 UTC (rev 2429)
+++ pkg/RcppDE/src/permuteMP.cpp 2010-11-12 03:28:41 UTC (rev 2430)
@@ -1,50 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
-//
-// Port of DEoptim (2.0.7) by Ardia et al to Rcpp/RcppArmadillo/Armadillo
-// Copyright (C) 2010 Dirk Eddelbuettel <edd at debian.org>
-//
-// DEoptim is Copyright (C) 2009 David Ardia and Katharine Mullen
-// and based on DE-Engine v4.0, Rainer Storn, 2004
-// (http://www.icsi.berkeley.edu/~storn/DeWin.zip)
-
-#ifdef USE_OPENMP
-#include <RcppArmadillo.h> // declarations for both Rcpp and RcppArmadillo offering Armadillo classes
-#include <omp.h> // OpenMP for compiler-generated multithreading
-
-// Function : void permute(int ia_urn2[], int i_urn2_depth)
-// Author : Rainer Storn (w/bug fixes contributed by DEoptim users)
-// Description : Generates i_urn2_depth random indices ex [0, i_NP-1]
-// which are all distinct. This is done by using a
-// permutation algorithm called the "urn algorithm"
-// which goes back to C.L.Robinson.
-// Functions : -
-// Globals : -
-// Parameters : ia_urn2 (O) array containing the random indices
-// i_urn2_depth (I) number of random indices (avoided index included)
-// i_NP (I) range of indices is [0, i_NP-1]
-// i_avoid (I) is the index to avoid and is located in ia_urn2[0].
-// ia_urn1 (I) additional temp vector
-// Preconditions : # Make sure that ia_urn2[] has a length of i_urn2_depth.
-// # i_urn2_depth must be smaller than i_NP.
-// Postconditions : # the index to be avoided is in ia_urn2[0], so fetch the
-// indices from ia_urn2[i], i = 1, 2, 3, ..., i_urn2_depth.
-// Return Value : -
-void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urn1[]) {
- //GetRNGstate();
- int i_urn1 = 0, i_urn2 = 0, k = i_NP;
- for (int i = 0; i < i_NP; i++)
- ia_urn1[i] = i; /* initialize urn1 */
- i_urn1 = i_avoid; /* get rid of the index to be avoided and place it in position 0. */
- // too simple pragma omp parallel for shared(ia_urn2,ia_urn1,i_urn1,i_urn2) private(k) schedule(static)
- // WORKS pragma omp parallel for shared(ia_urn2) private(k) schedule(dynamic)
- for (k = i_NP; k > i_NP - i_urn2_depth; k--) {
- ia_urn2[i_urn2] = ia_urn1[i_urn1]; /* move it into urn2 */
- ia_urn1[i_urn1] = ia_urn1[k-1]; /* move highest index to fill gap */
- //k = k - 1; /* reduce number of accessible indices */
- i_urn2 = i_urn2 + 1; /* next position in urn2 */
- //i_urn1 = static_cast<int>(::unif_rand() * k); /* choose a random index */
- i_urn1 = static_cast<int>(arma::randu() * (k-1)); /* choose a random index */
- }
- //PutRNGstate();
-}
-#endif
Copied: pkg/RcppParDE/DESCRIPTION (from rev 2425, pkg/RcppDE/DESCRIPTION)
===================================================================
--- pkg/RcppParDE/DESCRIPTION (rev 0)
+++ pkg/RcppParDE/DESCRIPTION 2010-11-12 03:28:41 UTC (rev 2430)
@@ -0,0 +1,17 @@
+Package: RcppParDE
+Version: 0.1.0
+Title: Global optimization by differential evolution in C++ using OpenMP parallel computing
+Author: Dirk Eddelbuettel extending DEoptim (by David Ardia, Katharine Mullen,
+ Brian Peterson, Joshua Ulrich) which itself is based on DE-Engine (by Rainer Storn)
+Maintainer: Dirk Eddelbuettel <edd at debian.org>
+Description: This package provides an efficient C++ based implementation of the
+ DEoptim function which performs global optimization by differential evolution.
+ It builds on RcppDE package which aims to show that "easier, shorter,
+ faster: pick any three" is achievable when moving code from plain old C to
+ modern C++ --- and attempts to reap further gains by using parallel
+ execution using the OpenMP framework for multithreaded computing on
+ multicore systems.
+License: GPL (>= 2)
+Depends: Rcpp, RcppArmadillo (>= 0.2.8)
+Suggests: inline, DEoptim
+LinkingTo: Rcpp, RcppArmadillo
Copied: pkg/RcppParDE/NAMESPACE (from rev 2425, pkg/RcppDE/NAMESPACE)
===================================================================
--- pkg/RcppParDE/NAMESPACE (rev 0)
+++ pkg/RcppParDE/NAMESPACE 2010-11-12 03:28:41 UTC (rev 2430)
@@ -0,0 +1,4 @@
+useDynLib(RcppParDE)
+export(DEoptim, DEoptim.control)
+S3method("plot", "DEoptim")
+S3method("summary", "DEoptim")
Copied: pkg/RcppParDE/R/DEoptim.R (from rev 2425, pkg/RcppDE/R/DEoptim.R)
===================================================================
--- pkg/RcppParDE/R/DEoptim.R (rev 0)
+++ pkg/RcppParDE/R/DEoptim.R 2010-11-12 03:28:41 UTC (rev 2430)
@@ -0,0 +1,152 @@
+DEoptim.control <- function(VTR = -Inf, strategy = 2, bs = FALSE, NP = 50,
+ itermax = 200, CR = 0.5, F = 0.8, trace = TRUE,
+ initialpop = NULL, storepopfrom = itermax + 1,
+ storepopfreq = 1, checkWinner = FALSE,
+ avWinner = TRUE, p = 0.2) {
+ if (itermax <= 0) {
+ warning("'itermax' <= 0; set to default value 200\n", immediate. = TRUE)
+ itermax <- 200
+ }
+ if (NP < 4) {
+ warning("'NP' < 4; set to default value 50\n", immediate. = TRUE)
+ NP <- 50
+ }
+ if (F < 0 | F > 2) {
+ warning("'F' not in [0,2]; set to default value 0.8\n", immediate. = TRUE)
+ F <- 0.8
+ }
+ if (CR < 0 | CR > 1) {
+ warning("'CR' not in [0,1]; set to default value 0.5\n", immediate. = TRUE)
+ CR <- 0.5
+ }
+ if (strategy < 1 | strategy > 6) {
+ warning("'strategy' not in {1,...,6}; set to default value 2\n",
+ immediate. = TRUE)
+ strategy <- 2
+ }
+
+ bs <- (bs > 0)
+
+ if ( trace < 0 ) {
+ warning("'trace' cannot be negative; set to 'TRUE'")
+ trace <- TRUE
+ }
+
+ storepopfreq <- floor(storepopfreq)
+ if (storepopfreq > itermax)
+ storepopfreq <- 1
+
+ if (p <= 0 || p > 1) {
+ warning("'p' not in (0,1]; set to default value 0.2\n", immediate. = TRUE)
+ p <- 0.2
+ }
+
+ list(VTR = VTR, strategy = strategy, NP = NP, itermax = itermax, CR
+ = CR, F = F, bs = bs, trace = trace, initialpop = initialpop,
+ storepopfrom = storepopfrom, storepopfreq = storepopfreq,
+ checkWinner = checkWinner, avWinner = avWinner, p = p)
+}
+
+DEoptim <- function(fn, lower, upper, control = DEoptim.control(), env, ...) {
+ ##fn1 <- function(par) fn(par, ...)
+ if (length(lower) != length(upper))
+ stop("'lower' and 'upper' are not of same length")
+ if (!is.vector(lower))
+ lower <- as.vector(lower)
+ if (!is.vector(upper))
+ upper <- as.vector(upper)
+ if (any(lower > upper))
+ stop("'lower' > 'upper'")
+ if (any(lower == "Inf"))
+ warning("you set a component of 'lower' to 'Inf'. May imply 'NaN' results", immediate. = TRUE)
+ if (any(lower == "-Inf"))
+ warning("you set a component of 'lower' to '-Inf'. May imply 'NaN' results", immediate. = TRUE)
+ if (any(upper == "Inf"))
+ warning("you set a component of 'upper' to 'Inf'. May imply 'NaN' results", immediate. = TRUE)
+ if (any(upper == "-Inf"))
+ warning("you set a component of 'upper' to '-Inf'. May imply 'NaN' results", immediate. = TRUE)
+ if (!is.null(names(lower)))
+ nam <- names(lower)
+ else if (!is.null(names(upper)) & is.null(names(lower)))
+ nam <- names(upper)
+ else
+ nam <- paste("par", 1:length(lower), sep = "")
+ if (missing(env))
+ env <- new.env()
+
+ ctrl <- do.call(DEoptim.control, as.list(control))
+ ctrl$npar <- length(lower)
+ if (ctrl$NP < 4) {
+ warning("'NP' < 4; set to default value 50\n", immediate. = TRUE)
+ ctrl$NP <- 50
+ }
+ if (ctrl$NP < 10*length(lower))
+ warning("For many problems it is best to set 'NP' (in 'control') to be at least ten times the length of the parameter vector. \n", immediate. = TRUE)
+ if (!is.null(ctrl$initialpop)) {
+ ctrl$specinitialpop <- TRUE
+ if(!identical(as.numeric(dim(ctrl$initialpop)), c(ctrl$NP, ctrl$npar)))
+ stop("Initial population is not a matrix with dim. NP x length(upper).")
+ }
+ else {
+ ctrl$specinitialpop <- FALSE
+ ctrl$initialpop <- matrix(0,1,1) # dummy matrix
+ }
+ ##
+ ctrl$trace <- as.numeric(ctrl$trace)
+ ctrl$specinitialpop <- as.numeric(ctrl$specinitialpop)
+
+ outC <- .Call("DEoptim", lower, upper, fn, ctrl, env, PACKAGE = "RcppParDE")
+ ##
+ if (length(outC$storepop) > 0) {
+ nstorepop <- floor((outC$iter - ctrl$storepopfrom) / ctrl$storepopfreq)
+ ## storepop <- list()
+ ## cnt <- 1
+ ## for(i in 1:nstorepop) {
+ ## idx <- cnt:((cnt - 1) + (ctrl$NP * ctrl$npar))
+ ## storepop[[i]] <- matrix(outC$storepop[idx], nrow = ctrl$NP, ncol = ctrl$npar,
+ ## byrow = TRUE)
+ ## cnt <- cnt + (ctrl$NP * ctrl$npar)
+ ## dimnames(storepop[[i]]) <- list(1:ctrl$NP, nam)
+ ## }
+ storepop <- outC$storepop[1:nstorepop]
+ for (i in 1:length(storepop)) dimnames(storepop[[i]]) <- list(1:ctrl$NP, nam)
+ }
+ else {
+ storepop = NULL
+ }
+
+ ## optim
+ bestmem <- as.numeric(outC$bestmem)
+ names(bestmem) <- nam
+ bestval <- as.numeric(outC$bestval)
+ nfeval <- as.numeric(outC$nfeval)
+ iter <- as.numeric(outC$iter)
+
+ ## member
+ names(lower) <- names(upper) <- nam
+ #bestmemit <- matrix(outC$bestmemit, nrow = iter, ncol = ctrl$npar, byrow = TRUE)
+ bestmemit <- outC$bestmemit
+
+ dimnames(bestmemit) <- list(1:iter, nam)
+ bestvalit <- as.numeric(outC$bestvalit[1:iter])
+ #pop <- matrix(outC$pop, nrow = ctrl$NP, ncol = ctrl$npar, byrow = TRUE)
+ pop <- outC$pop
+ storepop <- as.list(storepop)
+
+ outR <- list(optim = list(
+ bestmem = bestmem,
+ bestval = bestval,
+ nfeval = nfeval,
+ iter = iter),
+ member = list(
+ lower = lower,
+ upper = upper,
+ bestmemit = bestmemit,
+ bestvalit = bestvalit,
+ pop = pop,
+ storepop = storepop))
+
+ attr(outR, "class") <- "DEoptim"
+ return(outR)
+}
+
Copied: pkg/RcppParDE/R/methods.R (from rev 2425, pkg/RcppDE/R/methods.R)
===================================================================
--- pkg/RcppParDE/R/methods.R (rev 0)
+++ pkg/RcppParDE/R/methods.R 2010-11-12 03:28:41 UTC (rev 2430)
@@ -0,0 +1,64 @@
+'summary.DEoptim' <- function(object, ...){
+ digits <- max(5, getOption('digits') - 2)
+
+ cat("\n***** summary of DEoptim object *****",
+ "\nbest member : ", round(object$optim$bestmem, digits),
+ "\nbest value : ", round(object$optim$bestval, digits),
+ "\nafter : ", round(object$optim$iter), "generations",
+ "\nfn evaluated : ", round(object$optim$nfeval), "times",
+ "\n*************************************\n")
+
+ invisible(object)
+}
+
+plot.DEoptim <- function (x, plot.type = c("bestmemit", "bestvalit",
+"storepop"), ...) {
+ z <- x$member
+ niter <- length(z$bestvalit)
+ npar <- length(z$lower)
+ nam <- names(z$lower)
+ nstorepop <- length(z$storepop)
+ if (identical(plot.type[1], "bestmemit")) {
+ plot.new()
+ par(mfrow = c(min(3, npar), 1))
+ for (i in 1:npar) {
+ if (identical(i%%4, 0)) {
+ cat("new plot\n")
+ devAskNewPage(ask = TRUE)
+ }
+ plot(1:niter, z$bestmemit[, i], xlim = c(1, niter),
+ las = 1, xlab = "iteration", ylab = "value",
+ main = nam[i], ...)
+ abline(h = c(z$lower[i], z$upper[i]), col = "red")
+ }
+ }
+ else if (identical(plot.type[1], "bestvalit")) {
+ plot(1:niter, z$bestvalit, xlim = c(1, niter), las = 1,
+ xlab = "iteration", ylab = "function value", main =
+"convergence plot",
+ ...)
+ }
+ else if (identical(plot.type[1], "storepop") && nstorepop > 0) {
+ plot.new()
+ par(mfrow = c(min(3, npar), 1))
+ for (i in 1:npar) {
+ if (identical(i%%4, 0)) {
+ cat("new plot\n")
+ devAskNewPage(ask = TRUE)
+ }
+ tmp <- NULL
+ for (j in 1:nstorepop) {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 2430
More information about the Rcpp-commits
mailing list