[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