[Rcpp-commits] r2313 - pkg/RcppDE/src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sat Oct 16 21:21:31 CEST 2010


Author: edd
Date: 2010-10-16 21:21:27 +0200 (Sat, 16 Oct 2010)
New Revision: 2313

Modified:
   pkg/RcppDE/src/de4_0.cpp
   pkg/RcppDE/src/evaluate.cpp
Log:
no more pointer anywhere -- evaluate() now takes an arma row vector 
next is replacing all those tiny loops with row or column operations


Modified: pkg/RcppDE/src/de4_0.cpp
===================================================================
--- pkg/RcppDE/src/de4_0.cpp	2010-10-16 17:30:59 UTC (rev 2312)
+++ pkg/RcppDE/src/de4_0.cpp	2010-10-16 19:21:27 UTC (rev 2313)
@@ -2,48 +2,29 @@
 //
 // Port of DEoptim (2.0.7) 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)
 
-/*
-Implementation of DE based loosely on DE-Engine v4.0, Rainer Storn, 2004   
-by Katharine Mullen, 2009
-modified by Joshua Ulrich, 2010
-
-Storn's MS Visual C++ v5.0 was downloaded from 
-http://www.icsi.berkeley.edu/~storn/DeWin.zip
-
-Note that the struct t_pop was not used since we want to store in it a
-parameter vector of arbitrary length -- 
-
-using a construction like 
-typedef struct t_pop 
-{
-  double fa_cost[MAXCOST];    
-  double fa_vector[];         
-} t_pop;
-I have not been able to use Calloc or R_alloc to set aside memory for 
-type t_pop dynamically; I did not look into the bowels but believe this
-is likely a limitation of these functions.  
-*/
-
 #include <RcppArmadillo.h>
 
 RcppExport SEXP DEoptimC(SEXP lower, SEXP upper, SEXP fn, SEXP control, SEXP rho);
 void devol(double VTR, double f_weight, double fcross, int i_bs_flag, 
-           double *lower, double *upper, SEXP fcall, SEXP rho, int i_trace,
+           arma::colvec & lower, arma::colvec & upper, SEXP fcall, SEXP rho, int i_trace,
            int i_strategy, int i_D, int i_NP, int i_itermax,
-           double *initialpopv, int i_storepopfreq, int i_storepopfrom,
+           arma::colvec & initpopv, int i_storepopfreq, int i_storepopfrom,
            int i_specinitialpop, int i_check_winner, int i_av_winner,
            arma::mat    & ta_popP, arma::mat    & ta_oldP, arma::mat    & ta_newP, arma::rowvec & t_bestP, 
 	   arma::colvec & ta_popC, arma::colvec & ta_oldC, arma::colvec & ta_newC, double       & t_bestC,	
-           double *t_bestitP, double *t_tmpP, double *tempP,
-           arma::colvec & d_pop, arma::colvec & d_storepop, 
-	   arma::colvec & d_bestmemit, arma::colvec & d_bestvalit,
-           int *gi_iter, double i_pPct, long & l_nfeval);
+           arma::colvec & t_bestitP, arma::rowvec & t_tmpP, arma::colvec & tempP,
+           arma::colvec & d_pop, arma::colvec & d_storepop, arma::colvec & d_bestmemit, arma::colvec & d_bestvalit,
+           int & i_iterations, double i_pPct, long & l_nfeval);
 void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urntmp[]);
-RcppExport double evaluate(long &l_nfeval, double *param, SEXP par, SEXP fcall, SEXP env);
+RcppExport double evaluate(long &l_nfeval, const arma::rowvec & param, SEXP par, SEXP fcall, SEXP env);
 
 RcppExport SEXP DEoptimC(SEXP lowerS, SEXP upperS, SEXP fnS, SEXP controlS, SEXP rhoS) {
-    BEGIN_RCPP ;
+    BEGIN_RCPP ;	// macro to fill in try part of try/catch exception handler
 
     Rcpp::Function fn(fnS);						// function to mininise
     Rcpp::Environment rho(rhoS); 					// environment to do it in
@@ -68,6 +49,10 @@
     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 maxbound(f_upper.begin(), f_upper.size(), false);
+    arma::colvec initpopv(initialpopv.begin(), initialpopv.size(), false);
+
     arma::mat ta_popP(i_NP*2, i_D);    					// Data structures for parameter vectors 
     arma::mat ta_oldP(i_NP,   i_D);
     arma::mat ta_newP(i_NP,   i_D);
@@ -79,7 +64,7 @@
     double t_bestC; 
 
     arma::colvec t_bestitP(i_D);
-    arma::colvec t_tmpP(i_D); 
+    arma::rowvec t_tmpP(i_D); 
     arma::colvec tempP(i_D);
 
     int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
@@ -89,18 +74,13 @@
     arma::colvec d_bestvalit(i_itermax); 	 
     int i_iter = 0;
 
-    /*---optimization--------------------------------------*/
-    devol(VTR, f_weight, f_cross, i_bs_flag, f_lower.begin(), f_upper.begin(), Rcpp::wrap(fn), Rcpp::wrap(rho), i_trace,
-	  i_strategy, i_D, i_NP, i_itermax,
-	  initialpopv.begin(), 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.memptr(), t_tmpP.memptr(), tempP.memptr(),
-	  d_pop, d_storepop, d_bestmemit, d_bestvalit,
-	  &i_iter, i_pPct, l_nfeval);
-    /*---end optimization----------------------------------*/
+    // call actual Differential Evolution optimization given the parameters
+    devol(VTR, f_weight, f_cross, i_bs_flag, minbound, maxbound, Rcpp::wrap(fn), Rcpp::wrap(rho), i_trace,
+	  i_strategy, i_D, i_NP, i_itermax, initpopv, 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, tempP,
+	  d_pop, d_storepop, d_bestmemit, d_bestvalit, i_iter, i_pPct, l_nfeval);
 
+    // and return a named list to R
     return Rcpp::List::create(Rcpp::Named("bestmem")   = Rcpp::wrap(t_bestP),		// sexp_bestmem,
 			      Rcpp::Named("bestval")   = t_bestC,       // sexp_bestval,
 			      Rcpp::Named("nfeval")    = l_nfeval,   	// sexp_nfeval,
@@ -109,20 +89,21 @@
 			      Rcpp::Named("bestvalit") = Rcpp::wrap(d_bestvalit),	// sexp_bestvalit,
 			      Rcpp::Named("pop")       = Rcpp::wrap(d_pop),    		// sexp_pop,
 			      Rcpp::Named("storepop")  = Rcpp::wrap(d_storepop));   	// sexp_storepop)
-    END_RCPP
+    END_RCPP   // macro to fill in catch() part of try/catch exception handler
 }
 
+
 void devol(double VTR, double f_weight, double f_cross, int i_bs_flag,
-           double *fa_minbound, double *fa_maxbound, SEXP fcall, SEXP rho, int trace, 
+           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,
-           double *initialpopv, int i_storepopfrom, int i_storepopfreq, 
+           arma::colvec & initialpopv, 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::rowvec & t_bestP, arma::colvec & ta_popC, arma::colvec & ta_oldC, arma::colvec & ta_newC, 
 	   double & t_bestC,
-           double *t_bestitP, double *t_tmpP, double *tempP,
+           arma::colvec & t_bestitP, arma::rowvec & t_tmpP, arma::colvec & tempP,
            arma::colvec &d_pop, arma::colvec &d_storepop, arma::colvec & d_bestmemit, arma::colvec & d_bestvalit,
-           int *gi_iter, double i_pPct, long & l_nfeval) {
+           int & i_iterations, double i_pPct, long & l_nfeval) {
 
     const int urn_depth = 5;   			// 4 + one index to avoid 
     Rcpp::NumericVector par(i_D);  		// initialize parameter vector to pass to evaluate function 
@@ -131,17 +112,13 @@
     int ia_urn2[urn_depth];
     Rcpp::IntegerVector ia_urntmp(i_NP); 	// so that we don't need to re-allocated each time in permute
 
-    int i_nstorepop, i_xav;
-    i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
-  
-    int popcnt, bestacnt, same; /* lazy cnters */
-    
+    int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
+    int i_xav, popcnt, bestacnt, same; 		// lazy cnters 
     double f_jitter, f_dither, t_bestitC, t_tmpC, tmp_best; 
     
     arma::mat initialpop(i_NP, i_D); 
 
     int i_pbest;    				// vars for DE/current-to-p-best/1 
-
     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 
@@ -172,18 +149,15 @@
     }
     l_nfeval = 0;    				// number of function evaluations (this is an input via DEoptim.control, but we over-write it?) 
 
-    /*------Initialization-----------------------------*/
-    for (i = 0; i < i_NP; i++) {
-	if (i_specinitialpop <= 0) { 	// random initial member 
+    for (i = 0; i < i_NP; i++) {		// ------Initialization-----------------------------
+	if (i_specinitialpop <= 0) { 		// random initial member 
 	    for (j = 0; j < i_D; j++) {
 		ta_popP.at(i,j) = fa_minbound[j] + unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
 	    }
 	} else { /* or user-specified initial member */
 	    ta_popP.row(i) = initialpop.row(i);
 	} 
-	arma::rowvec r = ta_popP.row(i);
-	ta_popC[i] = evaluate(l_nfeval, r.memptr(), par, fcall, rho);
-
+	ta_popC[i] = evaluate(l_nfeval, ta_popP.row(i), par, fcall, rho);
 	if (i == 0 || ta_popC[i] <= t_bestC) {
 	    t_bestC = ta_popC[i];
 	    //for (j = 0; j < i_D; j++)  
@@ -192,12 +166,10 @@
 	}
     }
 
-    /*---assign pointers to current ("old") population---*/
-    ta_oldP = ta_popP;
+    ta_oldP = ta_popP;				// ---assign pointers to current ("old") population---
     ta_oldC = ta_popC;
   
-    /*------Iteration loop--------------------------------------------*/
-    int i_iter = 0;
+    int i_iter = 0;				// ------Iteration loop--------------------------------------------
     popcnt = 0;
     bestacnt = 0;
     i_xav = 1;
@@ -242,8 +214,9 @@
 	for (i = 0; i < i_NP; i++) {
 
 	    /*t_tmpP is the vector to mutate and eventually select*/
-	    for (j = 0; j < i_D; j++) 
-		t_tmpP[j] = ta_oldP.at(i,j);
+	    //for (j = 0; j < i_D; j++) 
+	    //	t_tmpP[j] = ta_oldP.at(i,j);
+	    t_tmpP = ta_oldP.row(i);
 	    t_tmpC = ta_oldC[i];
 
 	    permute(ia_urn2, urn_depth, i_NP, i, ia_urntmp.begin()); /* Pick 4 random and distinct */
@@ -381,9 +354,8 @@
 	    }
 
 	    /*------Trial mutation now in t_tmpP-----------------*/
-	    /* Evaluate mutant in t_tmpP[]*/
-	    t_tmpC = evaluate(l_nfeval, t_tmpP, par, fcall, rho); 
-	    
+	    t_tmpC = evaluate(l_nfeval, t_tmpP, par, fcall, rho); 	    // Evaluate mutant in t_tmpP[]
+
 	    /* note that i_bs_flag means that we will choose the
 	     *best NP vectors from the old and new population later*/
 	    if (t_tmpC <= ta_oldC[i] || i_bs_flag) {
@@ -466,9 +438,8 @@
 		}
 	    if(same && i_iter > 1)  {
 		i_xav++;
-		/* if re-evaluation of winner */
-		tmp_best = evaluate(l_nfeval, t_bestP.memptr(), par, fcall, rho);
-		
+		tmp_best = evaluate(l_nfeval, t_bestP, par, fcall, rho);			// if re-evaluation of winner 
+
 		/* possibly letting the winner be the average of all past generations */
 		if(i_av_winner)
 		    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);
@@ -483,8 +454,8 @@
 	    t_bestitP[j] = t_bestP[j];
 	t_bestitC = t_bestC;
 
-	if( trace > 0 ) {
-	    if( (i_iter % trace) == 0 ) {
+	if( i_trace > 0 ) {
+	    if( (i_iter % i_trace) == 0 ) {
 		Rprintf("Iteration: %d bestvalit: %f bestmemit:", i_iter, t_bestC);
 		for (j = 0; j < i_D; j++)
 		    Rprintf("%12.6f", t_bestP[j]);
@@ -501,7 +472,7 @@
 	    k++;
 	}
     }
-    *gi_iter = i_iter;
+    i_iterations = i_iter;
 
     PutRNGstate();
 }

Modified: pkg/RcppDE/src/evaluate.cpp
===================================================================
--- pkg/RcppDE/src/evaluate.cpp	2010-10-16 17:30:59 UTC (rev 2312)
+++ pkg/RcppDE/src/evaluate.cpp	2010-10-16 19:21:27 UTC (rev 2313)
@@ -1,25 +1,21 @@
 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
 //
-// Port of DEoptim (2.0.7) to Rcpp/RcppArmadillo/Armadillo
+// 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
 
-#include <Rcpp.h>
+#include <RcppArmadillo.h>
 
-/*------objective function---------------------------------------*/
-
-RcppExport double evaluate(long &l_nfeval, double *param, SEXP parS, SEXP fcall, SEXP env) {
-    Rcpp::NumericVector par(parS);
-    memcpy(par.begin(), param, par.size() * sizeof(double));
-
-    l_nfeval++;  		// increment function evaluation count 
-
-    SEXP fn = ::Rf_lang2(fcall, par); 			// this can be done with Rcpp 
+RcppExport double evaluate(long &l_nfeval, const arma::rowvec & param, SEXP parS, SEXP fcall, SEXP env) {
+    Rcpp::NumericVector par(parS); 			// access parS as numeric vector to fill it
+    std::copy(param.begin(), param.end(), par.begin()); // STL way of copying
+    SEXP fn = ::Rf_lang2(fcall, par); 			// this could be done with Rcpp 
     SEXP sexp_fvec = ::Rf_eval(fn, env);		// but is still a lot slower right now
     double f_result = Rcpp::as<double>(sexp_fvec);
-   
     if (ISNAN(f_result)) 
 	::Rf_error("NaN value of objective function! \nPerhaps adjust the bounds.");
-   
+    l_nfeval++;  					// increment function evaluation count 
    return(f_result); 
 }
 



More information about the Rcpp-commits mailing list