[Rcpp-commits] r2335 - pkg/RcppDE/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 18 21:15:58 CEST 2010
Author: edd
Date: 2010-10-18 21:15:57 +0200 (Mon, 18 Oct 2010)
New Revision: 2335
Modified:
pkg/RcppDE/src/de4_0.cpp
pkg/RcppDE/src/evaluate.cpp
Log:
rows and columns interchanged as columns can be accessed to much quicker
Modified: pkg/RcppDE/src/de4_0.cpp
===================================================================
--- pkg/RcppDE/src/de4_0.cpp 2010-10-18 19:01:14 UTC (rev 2334)
+++ pkg/RcppDE/src/de4_0.cpp 2010-10-18 19:15:57 UTC (rev 2335)
@@ -12,17 +12,15 @@
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,
- arma::rowvec & lower, arma::rowvec & 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,
arma::mat & initpopm, 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::rowvec & ta_popC, arma::rowvec & ta_oldC, arma::rowvec & ta_newC, double & t_bestC,
- arma::rowvec & t_bestitP, arma::rowvec & t_tmpP, arma::rowvec & tempP,
- arma::mat & d_pop, Rcpp::List & d_storepop, arma::mat & d_bestmemit, arma::rowvec & d_bestvalit,
+ 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::colvec & tempP,
+ 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);
-//RcppExport void permute(arma::icolvec & ia_urn2, int i_avoid, arma::icolvec & ia_urntmp);
-//RcppExport double evaluate(long & l_nfeval, const arma::rowvec & param, SEXP par, SEXP fcall, SEXP env);
void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urntmp[]);
RcppExport double evaluate(long &l_nfeval, const double *param, SEXP parS, SEXP fcall, SEXP env);
@@ -53,29 +51,29 @@
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::rowvec minbound(f_lower.begin(), f_lower.size(), false); // convert three Rcpp vectors to arma vectors
- arma::rowvec maxbound(f_upper.begin(), f_upper.size(), false);
+ 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::mat initpopm(initialpopm.begin(), initialpopm.rows(), initialpopm.cols(), 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);
- arma::rowvec t_bestP(i_D);
+ 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::rowvec ta_popC(i_NP*2); // Data structures for obj. fun. values associated with par. vectors
- arma::rowvec ta_oldC(i_NP);
- arma::rowvec ta_newC(i_NP);
+ arma::colvec ta_popC(i_NP*2); // Data structures for obj. fun. values associated with par. vectors
+ arma::colvec ta_oldC(i_NP);
+ arma::colvec ta_newC(i_NP);
double t_bestC;
- arma::rowvec t_bestitP(i_D);
- arma::rowvec t_tmpP(i_D);
- arma::rowvec tempP(i_D);
+ arma::colvec t_bestitP(i_D);
+ arma::colvec t_tmpP(i_D);
+ arma::colvec tempP(i_D);
int i_nstorepop = ceil((i_itermax - i_storepopfrom) / i_storepopfreq);
- arma::mat d_pop(i_NP, i_D);
+ arma::mat d_pop(i_D, i_NP);
Rcpp::List d_storepop(i_nstorepop);
- arma::mat d_bestmemit(i_itermax, i_D);
- arma::rowvec d_bestvalit(i_itermax);
+ arma::mat d_bestmemit(i_D, i_itermax);
+ arma::colvec d_bestvalit(i_itermax);
int i_iter = 0;
// call actual Differential Evolution optimization given the parameters
@@ -89,50 +87,44 @@
Rcpp::Named("bestval") = t_bestC, // sexp_bestval,
Rcpp::Named("nfeval") = l_nfeval, // sexp_nfeval,
Rcpp::Named("iter") = i_iter, // sexp_iter,
- Rcpp::Named("bestmemit") = d_bestmemit, // sexp_bestmemit,
+ Rcpp::Named("bestmemit") = trans(d_bestmemit), // sexp_bestmemit,
Rcpp::Named("bestvalit") = d_bestvalit, // sexp_bestvalit,
- Rcpp::Named("pop") = d_pop, // sexp_pop,
+ Rcpp::Named("pop") = trans(d_pop), // sexp_pop,
Rcpp::Named("storepop") = d_storepop); // sexp_storepop)
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,
- arma::rowvec & fa_minbound, arma::rowvec & fa_maxbound, SEXP fcall, SEXP rho, int i_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, 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::rowvec & t_bestP,
- arma::rowvec & ta_popC, arma::rowvec & ta_oldC, arma::rowvec & ta_newC, double & t_bestC,
- arma::rowvec & t_bestitP, arma::rowvec & t_tmpP, arma::rowvec & tempP,
- arma::mat &d_pop, Rcpp::List &d_storepop, arma::mat & d_bestmemit, arma::rowvec & d_bestvalit,
+ 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::colvec & tempP,
+ 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) {
const int urn_depth = 5; // 4 + one index to avoid
- //arma::rowvec parvec(i_D); // initialize parameter vector to pass to evaluate function
- //SEXP par = Rcpp::wrap(parvec);
- Rcpp::NumericVector par(i_D);
+ Rcpp::NumericVector par(i_D); // initialize parameter vector to pass to evaluate function
int i, j, k, i_r1, i_r2, i_r3, i_r4; // counting variables and placeholders for random indexes
- int ia_urn2[urn_depth];
- //Rcpp::IntegerVector ia_urntmp(i_NP); // so that we don't need to re-allocated each time in permute
- arma::irowvec ia_urntmp(i_NP);
- //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::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
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;
+ double f_jitter, f_dither, t_bestitC, t_tmpC, tmp_best; // , tempC
- arma::mat initialpop(i_NP, i_D);
+ arma::mat initialpop(i_D, i_NP);
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::irowvec sortIndex(i_NP); // sorted values of ta_oldC
+ arma::icolvec sortIndex(i_NP); // sorted values of ta_oldC
for(i = 0; i < i_NP; i++) sortIndex[i] = i;
int i_len, done, step, bound; // vars for when i_bs_flag == 1 */
- double tempC;
GetRNGstate();
@@ -143,23 +135,22 @@
i_nstorepop = (i_nstorepop < 0) ? 0 : i_nstorepop;
if (i_specinitialpop > 0) { // if initial population provided, initialize with values
- initialpop = initialpopm;
+ initialpop = trans(initialpopm); // transpose as we prefer columns for population members here
}
//l_nfeval = 0; // already init'ed in main function
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]);
+ ta_popP.at(j,i) = fa_minbound[j] + ::unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
}
} else { // or user-specified initial member
- ta_popP.row(i) = initialpop.row(i);
+ ta_popP.col(i) = initialpop.col(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.colptr(i), par, fcall, rho);
if (i == 0 || ta_popC[i] <= t_bestC) {
t_bestC = ta_popC[i];
- t_bestP = ta_popP.row(i);
+ t_bestP = ta_popP.unsafe_col(i);
}
}
@@ -173,10 +164,10 @@
while ((i_iter < i_itermax) && (t_bestC > VTR)) { // main loop ====================================
if (i_iter % i_storepopfreq == 0 && i_iter >= i_storepopfrom) { // store intermediate populations
- d_storepop[popcnt++] = Rcpp::wrap(ta_oldP);
+ d_storepop[popcnt++] = Rcpp::wrap( trans(ta_oldP) );
} // end store pop
- d_bestmemit.row(i_iter) = t_bestP; // store the best member
+ 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;
t_bestitC = t_bestC;
@@ -185,17 +176,17 @@
f_dither = f_weight + ::unif_rand() * (1.0 - f_weight); // ----computer dithering factor -----------------
if (i_strategy == 6) { // ---DE/current-to-p-best/1 -----------------------------------------------------
- arma::rowvec temp_oldC = ta_oldC; // create a copy of ta_oldC to avoid changing it
+ arma::colvec temp_oldC = ta_oldC; // create a copy of ta_oldC to avoid changing it
rsort_with_index( temp_oldC.memptr(), sortIndex.begin(), i_NP ); // sort temp_oldC to use sortIndex later
}
for (i = 0; i < i_NP; i++) { // ----start of loop through ensemble------------------------
- t_tmpP = ta_oldP.row(i); // t_tmpP is the vector to mutate and eventually select
+ t_tmpP = ta_oldP.col(i); // t_tmpP is the vector to mutate and eventually select
t_tmpC = ta_oldC[i];
- permute(ia_urn2, urn_depth, i_NP, i, ia_urntmp.memptr()); // Pick 4 random and distinct
//permute(ia_urn2, i, ia_urntmp); // Pick 4 random and distinct
+ permute(ia_urn2.memptr(), urn_depth, i_NP, i, ia_urntmp.memptr()); // Pick 4 random and distinct
i_r1 = ia_urn2[1]; // population members
i_r2 = ia_urn2[2];
i_r3 = ia_urn2[3];
@@ -208,7 +199,7 @@
case 2: // ---DE/local-to-best/1/bin---------------------------------------------------
j = static_cast<int>(::unif_rand() * 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(i_r2,j) - ta_oldP.at(i_r3,j));
+ t_tmpP[j] = t_tmpP[j] + f_weight * (t_bestitP[j] - t_tmpP[j]) + f_weight * (ta_oldP.at(j,i_r2) - ta_oldP.at(j,i_r3));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -216,7 +207,7 @@
case 1: // ---classical strategy DE/rand/1/bin-----------------------------------------
j = static_cast<int>(::unif_rand() * i_D); // random parameter
do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(i_r1,j) + f_weight * (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
+ t_tmpP[j] = ta_oldP.at(j,i_r1) + f_weight * (ta_oldP.at(j,i_r2) - ta_oldP.at(j,i_r3));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -225,7 +216,7 @@
j = static_cast<int>(::unif_rand() * i_D); // random parameter
do { // add fluctuation to random target
f_jitter = 0.0001 * ::unif_rand() + f_weight;
- t_tmpP[j] = t_bestitP[j] + f_jitter * (ta_oldP.at(i_r1,j) - ta_oldP.at(i_r2,j));
+ t_tmpP[j] = t_bestitP[j] + f_jitter * (ta_oldP.at(j,i_r1) - ta_oldP.at(j,i_r2));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -233,7 +224,7 @@
case 4: // ---DE/rand/1/bin with per-vector-dither-------------------------------------
j = static_cast<int>(::unif_rand() * i_D); // random parameter
do { // add fluctuation to random target *
- t_tmpP[j] = ta_oldP.at(i_r1,j) + (f_weight + ::unif_rand()*(1.0 - f_weight))* (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
+ t_tmpP[j] = ta_oldP.at(j,i_r1) + (f_weight + ::unif_rand()*(1.0 - f_weight))* (ta_oldP.at(j,i_r2) - ta_oldP.at(j,i_r3));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -241,7 +232,7 @@
case 5: // ---DE/rand/1/bin with per-generation-dither---------------------------------
j = static_cast<int>(::unif_rand() * i_D); // random parameter
do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(i_r1,j) + f_dither * (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
+ t_tmpP[j] = ta_oldP.at(j,i_r1) + f_dither * (ta_oldP.at(j,i_r2) - ta_oldP.at(j,i_r3));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -250,7 +241,7 @@
i_pbest = sortIndex[static_cast<int>(::unif_rand() * p_NP)]; // select from [0, 1, 2, ..., (pNP-1)]
j = static_cast<int>(::unif_rand() * i_D); // random parameter
do { // add fluctuation to random target
- t_tmpP[j] = ta_oldP.at(i,j) + f_weight * (ta_oldP.at(i_pbest,j) - ta_oldP.at(i,j)) + f_weight * (ta_oldP.at(i_r1,j) - ta_oldP.at(i_r2,j));
+ 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,i_r1) - ta_oldP.at(j,i_r2));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
break;
@@ -258,16 +249,14 @@
default: // ---variation to DE/rand/1/bin: either-or-algorithm--------------------------
j = static_cast<int>(::unif_rand() * i_D); // random parameter
if (::unif_rand() < 0.5) { // differential mutation, Pmu = 0.5
- do {
- // add fluctuation to random target */
- t_tmpP[j] = ta_oldP.at(i_r1,j) + f_weight * (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
+ do { // add fluctuation to random target */
+ t_tmpP[j] = ta_oldP.at(j,i_r1) + f_weight * (ta_oldP.at(j,i_r2) - ta_oldP.at(j,i_r3));
j = (j + 1) % i_D;
} while ((::unif_rand() < 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(i_r1,j) + 0.5 * (f_weight + 1.0) * (ta_oldP.at(i_r2,j) + ta_oldP.at(i_r3,j) - 2 * ta_oldP.at(i_r1,j));
+ do { // add fluctuation to random target */
+ t_tmpP[j] = ta_oldP.at(j,i_r1) + 0.5 * (f_weight + 1.0) * (ta_oldP.at(j,i_r2) + ta_oldP.at(j,i_r3) - 2 * ta_oldP.at(j,i_r1));
j = (j + 1) % i_D;
} while ((::unif_rand() < f_cross) && (++k < i_D));
}
@@ -286,16 +275,15 @@
// ------Trial mutation now in t_tmpP-----------------
t_tmpC = evaluate(l_nfeval, t_tmpP.memptr(), 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) {
- ta_newP.row(i) = t_tmpP; // replace target with mutant
+ if (t_tmpC <= ta_oldC[i] || i_bs_flag) { // i_bs_flag means that we will choose best NP vectors from old and new population 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.row(i) = ta_oldP.row(i);
+ ta_newP.col(i) = ta_oldP.col(i);
ta_newC[i] = ta_oldC[i];
}
} // End mutation loop through pop., ie the "for (i = 0; i < i_NP; i++)"
@@ -318,12 +306,8 @@
for (j = 0; j < bound; j++) {
i = j + step + 1;
if (ta_popC[j] > ta_popC[i-1]) {
- tempP = ta_popP.row(i-1);
- tempC = ta_popC[i-1];
- ta_popP.row(i-1) = ta_popP.row(j);
- ta_popC[i-1] = ta_popC[j];
- ta_popP.row(j) = tempP;
- ta_popC[j] = tempC;
+ ta_popP.swap_cols(i-1, j);
+ ta_popC.swap_rows(i-1, j);
done = 0;
} // if
} // for
@@ -373,7 +357,6 @@
}
inline void permute(int ia_urn2[], int i_urn2_depth, int i_NP, int i_avoid, int ia_urn1[])
-//RcppExport inline void permute(arma::icolvec & ia_urn2, int i_avoid, arma::icolvec & ia_urn1)
/********************************************************************
** Function : void permute(int ia_urn2[], int i_urn2_depth)
** Author : Rainer Storn (w/bug fixes contributed by DEoptim users)
@@ -397,9 +380,6 @@
{
GetRNGstate();
int k = i_NP;
- //int i_urn2_depth = ia_urn2.n_elem; // ie const urn_depth
- //int i_NP = ia_urn1.n_elem;
- //int k = ia_urn1.n_elem;
int i_urn1 = 0;
int i_urn2 = 0;
for (int i = 0; i < i_NP; i++)
Modified: pkg/RcppDE/src/evaluate.cpp
===================================================================
--- pkg/RcppDE/src/evaluate.cpp 2010-10-18 19:01:14 UTC (rev 2334)
+++ pkg/RcppDE/src/evaluate.cpp 2010-10-18 19:15:57 UTC (rev 2335)
@@ -10,7 +10,6 @@
//RcppExport double evaluate(long &l_nfeval, const arma::rowvec & param, SEXP parS, SEXP fcall, SEXP env) {
RcppExport double evaluate(long &l_nfeval, const double *param, SEXP parS, SEXP fcall, SEXP env) {
Rcpp::NumericVector par(parS); // access parS as numeric vector to fill it
- //memcpy(par.begin(), param.memptr(), par.size() * sizeof(double));
//std::copy(param.begin(), param.end(), par.begin()); // STL way of copying
memcpy(par.begin(), param, par.size() * sizeof(double));
SEXP fn = ::Rf_lang2(fcall, par); // this could be done with Rcpp
More information about the Rcpp-commits
mailing list