[Rcpp-commits] r2312 - in pkg/RcppDE: . src

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


Author: edd
Date: 2010-10-16 19:30:59 +0200 (Sat, 16 Oct 2010)
New Revision: 2312

Modified:
   pkg/RcppDE/benchmark.txt
   pkg/RcppDE/src/de4_0.cpp
Log:
t_bestP is now a rowvec so that we can start to copy row-wise at once
another re-indentation


Modified: pkg/RcppDE/benchmark.txt
===================================================================
--- pkg/RcppDE/benchmark.txt	2010-10-16 16:50:44 UTC (rev 2311)
+++ pkg/RcppDE/benchmark.txt	2010-10-16 17:30:59 UTC (rev 2312)
@@ -79,3 +79,19 @@
 Genrose20   0.887667 0.879889          0.99124      0.876205
 Genrose50   2.523611 2.518944          0.99815      0.184920
 MEANS       0.563456 0.561589          0.99669      0.331289
+
+
+# At  2010-10-16 12:26:24.650366 
+# SVN  2311M 
+             DEoptim   RcppDE ratioRcppToBasic pctGainOfRcpp
+Rastrigin2  0.040389 0.041333          1.02338      -2.33838
+Rastrigin5  0.105833 0.107833          1.01890      -1.88976
+Rastrigin20 0.543667 0.547778          1.00756      -0.75618
+Wild2       0.066667 0.066667          1.00000       0.00000
+Wild5       0.181500 0.181222          0.99847       0.15305
+Wild20      1.042833 1.039667          0.99696       0.30366
+Genrose2    0.070667 0.071778          1.01572      -1.57233
+Genrose5    0.186167 0.187222          1.00567      -0.56699
+Genrose20   0.894167 0.896722          1.00286      -0.28580
+Genrose50   2.527278 2.562778          1.01405      -1.40467
+MEANS       0.565917 0.570300          1.00775      -0.77455

Modified: pkg/RcppDE/src/de4_0.cpp
===================================================================
--- pkg/RcppDE/src/de4_0.cpp	2010-10-16 16:50:44 UTC (rev 2311)
+++ pkg/RcppDE/src/de4_0.cpp	2010-10-16 17:30:59 UTC (rev 2312)
@@ -33,7 +33,7 @@
            int i_strategy, int i_D, int i_NP, int i_itermax,
            double *initialpopv, 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::colvec & t_bestP, 
+           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, 
@@ -71,7 +71,7 @@
     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::colvec t_bestP(i_D); 						// double *t_bestP = (double *)R_alloc(1,sizeof(double) * i_D);
+    arma::rowvec t_bestP(i_D); 
 
     arma::colvec ta_popC(i_NP*2);  				    	// Data structures for obj. fun. values associated with par. vectors 
     arma::colvec ta_oldC(i_NP);
@@ -118,7 +118,7 @@
            double *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::colvec & t_bestP, arma::colvec & ta_popC, arma::colvec & ta_oldC, arma::colvec & ta_newC, 
+	   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,
@@ -174,11 +174,12 @@
 
     /*------Initialization-----------------------------*/
     for (i = 0; i < i_NP; i++) {
-	for (j = 0; j < i_D; j++) {
-	    if (i_specinitialpop <= 0) { 	// random initial member 
+	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.at(i,j) = initialpop.at(i,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);
@@ -201,303 +202,297 @@
     bestacnt = 0;
     i_xav = 1;
   
-    /* loop */
-    while ((i_iter < i_itermax) && (t_bestC > VTR))
-	{
-	    /* store intermediate populations */
-	    if (i_iter % i_storepopfreq == 0 && i_iter >= i_storepopfrom) {
-		for (i = 0; i < i_NP; i++) {
-		    for (j = 0; j < i_D; j++) {
-			d_storepop[popcnt] = ta_oldP.at(i,j);
-			popcnt++;
-		    }
+    while ((i_iter < i_itermax) && (t_bestC > VTR)) {    // loop 
+	/* store intermediate populations */
+	if (i_iter % i_storepopfreq == 0 && i_iter >= i_storepopfrom) {
+	    for (i = 0; i < i_NP; i++) {
+		for (j = 0; j < i_D; j++) {
+		    d_storepop[popcnt] = ta_oldP.at(i,j);
+		    popcnt++;
 		}
-	    } /* end store pop */
-      
-	    /* store the best member */
-	    for(j = 0; j < i_D; j++) {
-		d_bestmemit[bestacnt] = t_bestP[j];
-		bestacnt++;
 	    }
-	    /* store the best value */
-	    d_bestvalit[i_iter] = t_bestC;
+	} /* end store pop */
       
-	    for (j = 0; j < i_D; j++) 
-		t_bestitP[j] = t_bestP[j];
-	    t_bestitC = t_bestC;
+	/* store the best member */
+	for(j = 0; j < i_D; j++) {
+	    d_bestmemit[bestacnt] = t_bestP[j];
+	    bestacnt++;
+	}
+	/* store the best value */
+	d_bestvalit[i_iter] = t_bestC;
+	
+	for (j = 0; j < i_D; j++) 
+	    t_bestitP[j] = t_bestP[j];
+	t_bestitC = t_bestC;
       
-	    i_iter++;
+	i_iter++;
      
-	    /*----computer dithering factor -----------------*/
-	    f_dither = f_weight + unif_rand() * (1.0 - f_weight);
+	/*----computer dithering factor -----------------*/
+	f_dither = f_weight + unif_rand() * (1.0 - f_weight);
       
-	    /*---DE/current-to-p-best/1 -----------------------------------------------------*/
-	    if (i_strategy == 6) {
-		/* create a copy of ta_oldC to avoid changing it */
-		arma::colvec temp_oldC = ta_oldC;
-		/* sort temp_oldC to use sortIndex later */
-		rsort_with_index( temp_oldC.memptr(), sortIndex.begin(), i_NP );
-	    }
+	/*---DE/current-to-p-best/1 -----------------------------------------------------*/
+	if (i_strategy == 6) {
+	    /* create a copy of ta_oldC to avoid changing it */
+	    arma::colvec temp_oldC = ta_oldC;
+	    /* sort temp_oldC to use sortIndex later */
+	    rsort_with_index( temp_oldC.memptr(), sortIndex.begin(), i_NP );
+	}
 
-	    /*----start of loop through ensemble------------------------*/
-	    for (i = 0; i < i_NP; i++) {
+	/*----start of loop through ensemble------------------------*/
+	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);
-		t_tmpC = ta_oldC[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);
+	    t_tmpC = ta_oldC[i];
 
-		permute(ia_urn2, urn_depth, i_NP, i, ia_urntmp.begin()); /* Pick 4 random and distinct */
-
-		i_r1 = ia_urn2[1];  /* population members */
-		i_r2 = ia_urn2[2];
-		i_r3 = ia_urn2[3];
-		i_r4 = ia_urn2[4];
-	
-		/*===Choice of strategy=======================================================*/
-		/*---classical strategy DE/rand/1/bin-----------------------------------------*/
-		if (i_strategy == 1) {
+	    permute(ia_urn2, urn_depth, i_NP, i, ia_urntmp.begin()); /* Pick 4 random and distinct */
+	    i_r1 = ia_urn2[1];  /* population members */
+	    i_r2 = ia_urn2[2];
+	    i_r3 = ia_urn2[3];
+	    i_r4 = ia_urn2[4];
+		
+	    /*===Choice of strategy=======================================================*/
+	    /*---classical strategy DE/rand/1/bin-----------------------------------------*/
+	    if (i_strategy == 1) {
 	  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
-		    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));
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
+		    
+		    j = (j + 1) % i_D;
+		    k++;
+		}while((unif_rand() < f_cross) && (k < i_D));
 
-			j = (j + 1) % i_D;
-			k++;
-		    }while((unif_rand() < f_cross) && (k < i_D));
-
-		}
-		/*---DE/local-to-best/1/bin---------------------------------------------------*/
-		else if (i_strategy == 2) {
+	    }
+	    /*---DE/local-to-best/1/bin---------------------------------------------------*/
+	    else if (i_strategy == 2) {
 	 
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
-		    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));
-			j = (j + 1) % i_D;
-			k++;
-		    }while((unif_rand() < f_cross) && (k < i_D));
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
+		    j = (j + 1) % i_D;
+		    k++;
+		} while((unif_rand() < f_cross) && (k < i_D));
 
-		}
-		/*---DE/best/1/bin with jitter------------------------------------------------*/
-		else if (i_strategy == 3) {
+	    }
+	    /*---DE/best/1/bin with jitter------------------------------------------------*/
+	    else if (i_strategy == 3) {
 	 	  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
-		    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));
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
 	    
-			j = (j + 1) % i_D;
-			k++;
-		    }while((unif_rand() < f_cross) && (k < i_D));
+		    j = (j + 1) % i_D;
+		    k++;
+		} while((unif_rand() < f_cross) && (k < i_D));
 
-		}
-		/*---DE/rand/1/bin with per-vector-dither-------------------------------------*/
-		else if (i_strategy == 4) {
+	    }
+	    /*---DE/rand/1/bin with per-vector-dither-------------------------------------*/
+	    else if (i_strategy == 4) {
 		  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
-		    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));
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
+		    
+		    j = (j + 1) % i_D;
+		    k++;
+		} while((unif_rand() < f_cross) && (k < i_D));
 
-			j = (j + 1) % i_D;
-			k++;
-		    }while((unif_rand() < f_cross) && (k < i_D));
+	    }
+	    /*---DE/rand/1/bin with per-generation-dither---------------------------------*/
+	    else if (i_strategy == 5) {
+	  
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
+	    
+		    j = (j + 1) % i_D;
+		    k++;
+		} while((unif_rand() < f_cross) && (k < i_D));
+       
+	    }
+	    /*---DE/current-to-p-best/1 (JADE)--------------------------------------------*/
+	    else if (i_strategy == 6) {
 
-		}
-		/*---DE/rand/1/bin with per-generation-dither---------------------------------*/
-		else if (i_strategy == 5) {
+		/* select from [0, 1, 2, ..., (pNP-1)] */
+		i_pbest = sortIndex[(int)(unif_rand() * p_NP)];
+		
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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));
+	    
+		    j = (j + 1) % i_D;
+		    k++;
+		}while((unif_rand() < f_cross) && (k < i_D));
+
+	    }
+	    /*---variation to DE/rand/1/bin: either-or-algorithm--------------------------*/
+	    else {
 	  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
+		j = (int)(unif_rand() * i_D); /* random parameter */
+		k = 0;
+		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_dither * (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
-	    
+			t_tmpP[j] = ta_oldP.at(i_r1,j) + f_weight * (ta_oldP.at(i_r2,j) - ta_oldP.at(i_r3,j));
+			
 			j = (j + 1) % i_D;
 			k++;
 		    }while((unif_rand() < f_cross) && (k < i_D));
-       
 		}
-		/*---DE/current-to-p-best/1 (JADE)--------------------------------------------*/
-		else if (i_strategy == 6) {
-
-		    /* select from [0, 1, 2, ..., (pNP-1)] */
-		    i_pbest = sortIndex[(int)(unif_rand() * p_NP)];
-	  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
+		else {
+		    /* recombination with K = 0.5*(F+1) -. F-K-Rule */
 		    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(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));
+			
 			j = (j + 1) % i_D;
 			k++;
 		    }while((unif_rand() < f_cross) && (k < i_D));
 
 		}
-		/*---variation to DE/rand/1/bin: either-or-algorithm--------------------------*/
-		else {
-	  
-		    j = (int)(unif_rand() * i_D); /* random parameter */
-		    k = 0;
-		    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));
-	      
-			    j = (j + 1) % i_D;
-			    k++;
-			}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));
-	      
-			    j = (j + 1) % i_D;
-			    k++;
-			}while((unif_rand() < f_cross) && (k < i_D));
-
-		    }
-		}/* end if (i_strategy ...*/
+	    }/* end if (i_strategy ...*/
 	
-		/*----boundary constraints, bounce-back method was not enforcing bounds correctly*/
-		for (j = 0; j < i_D; j++) {
-		    if (t_tmpP[j] < fa_minbound[j]) {
-			t_tmpP[j] = fa_minbound[j] +
-			    unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
-		    }
-		    if (t_tmpP[j] > fa_maxbound[j]) {
-			t_tmpP[j] =  fa_maxbound[j] -
-			    unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
-		    }
+	    /*----boundary constraints, bounce-back method was not enforcing bounds correctly*/
+	    for (j = 0; j < i_D; j++) {
+		if (t_tmpP[j] < fa_minbound[j]) {
+		    t_tmpP[j] = fa_minbound[j] +
+			unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
 		}
+		if (t_tmpP[j] > fa_maxbound[j]) {
+		    t_tmpP[j] =  fa_maxbound[j] -
+			unif_rand() * (fa_maxbound[j] - fa_minbound[j]);
+		}
+	    }
 
-		/*------Trial mutation now in t_tmpP-----------------*/
-		/* Evaluate mutant in t_tmpP[]*/
-		t_tmpC = evaluate(l_nfeval, t_tmpP, par, fcall, rho); 
-	
-		/* 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) {
-		    /* replace target with mutant */
+	    /*------Trial mutation now in t_tmpP-----------------*/
+	    /* Evaluate mutant in t_tmpP[]*/
+	    t_tmpC = evaluate(l_nfeval, t_tmpP, par, fcall, rho); 
+	    
+	    /* 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) {
+		/* replace target with mutant */
+		for (j = 0; j < i_D; j++) 
+		    ta_newP.at(i,j) = t_tmpP[j];
+		ta_newC[i] = t_tmpC;
+		if (t_tmpC <= t_bestC) {
 		    for (j = 0; j < i_D; j++) 
-			ta_newP.at(i,j) = t_tmpP[j];
-		    ta_newC[i] = t_tmpC;
-		    if (t_tmpC <= t_bestC) {
-			for (j = 0; j < i_D; j++) 
-			    t_bestP[j] = t_tmpP[j];
-			t_bestC = t_tmpC;
-		    }
-		} 
-		else {
-		    for (j = 0; j < i_D; j++) 
-			ta_newP.at(i,j) = ta_oldP.at(i,j);
-		    ta_newC[i] = ta_oldC[i];
-	  
+			t_bestP[j] = t_tmpP[j];
+		    t_bestC = t_tmpC;
 		}
-	    } /* End mutation loop through pop. */
+	    } 
+	    else {
+		for (j = 0; j < i_D; j++) 
+		    ta_newP.at(i,j) = ta_oldP.at(i,j);
+		ta_newC[i] = ta_oldC[i];
+	  
+	    }
+	} /* End mutation loop through pop. */
+	
      
-     
-	    if(i_bs_flag) {
-		/* examine old and new pop. and take the best NP members
-		 * into next generation */
-		for (i = 0; i < i_NP; i++) {
-		    for (j = 0; j < i_D; j++) 
-			ta_popP.at(i,j) = ta_oldP.at(i,j);
-		    ta_popC[i] = ta_oldC[i];
-		}
-		for (i = 0; i < i_NP; i++) {
-		    for (j = 0; j < i_D; j++) 
-			ta_popP.at(i_NP+i,j) = ta_newP.at(i,j);
-		    ta_popC[i_NP+i] = ta_newC[i];
-		}
-		i_len = 2 * i_NP;
-		step = i_len;  /* array length */
-		while (step > 1) {
-		    step /= 2;   /* halve the step size */
-		    do {
-			done = 1;
-			bound  = i_len - step;
-			for (j = 0; j < bound; j++) {
-			    i = j + step + 1;
-			    if (ta_popC[j] > ta_popC[i-1]) {
-				for (k = 0; k < i_D; k++) 
-				    tempP[k] = ta_popP.at(i-1, k);
-				tempC = ta_popC[i-1];
-				for (k = 0; k < i_D; k++) 
-				    ta_popP.at(i-1,k) = ta_popP.at(j,k);
-				ta_popC[i-1] = ta_popC[j];
-				for (k = 0; k < i_D; k++) 
-				    ta_popP.at(j,k) = tempP[k];
-				ta_popC[j] = tempC;
-				done = 0; 
-				/* if a swap has been made we are not finished yet */
-			    }  /* if */
-			}  /* for */
-		    } while (!done);   /* while */
-		} /*while (step > 1) */
-		/* now the best NP are in first NP places in gta_pop, use them */
-		for (i = 0; i < i_NP; i++) {
-		    for (j = 0; j < i_D; j++) 
-			ta_newP.at(i,j) = ta_popP.at(i,j);
-		    ta_newC[i] = ta_popC[i];
-		}
-	    } /*i_bs_flag*/
-
-	    /* have selected NP mutants move on to next generation */
+	if(i_bs_flag) {
+	    /* examine old and new pop. and take the best NP members into next generation */
 	    for (i = 0; i < i_NP; i++) {
 		for (j = 0; j < i_D; j++) 
-		    ta_oldP.at(i,j) = ta_newP.at(i,j);
-		ta_oldC[i] = ta_newC[i];
+		    ta_popP.at(i,j) = ta_oldP.at(i,j);
+		ta_popC[i] = ta_oldC[i];
 	    }
-	    /* check if the best stayed the same, if necessary */
-	    if(i_check_winner)  {
-		same = 1;
-		for (j = 0; j < i_D; j++)
-		    if (t_bestitP[j] != t_bestP[j]) {
-			same = 0;
-		    }
-		if(same && i_iter > 1)  {
-		    i_xav++;
-		    /* if re-evaluation of winner */
-		    tmp_best = evaluate(l_nfeval, t_bestP.memptr(), par, fcall, rho);
-	 
-		    /* 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);
-		    else
-			t_bestC = tmp_best;
-	
-		}
-		else {
-		    i_xav = 1;
-		}
-	
+	    for (i = 0; i < i_NP; i++) {
+		for (j = 0; j < i_D; j++) 
+		    ta_popP.at(i_NP+i,j) = ta_newP.at(i,j);
+		ta_popC[i_NP+i] = ta_newC[i];
 	    }
-	    for (j = 0; j < i_D; j++) 
-		t_bestitP[j] = t_bestP[j];
-	    t_bestitC = t_bestC;
+	    i_len = 2 * i_NP;
+	    step = i_len;  /* array length */
+	    while (step > 1) {
+		step /= 2;   /* halve the step size */
+		do {
+		    done = 1;
+		    bound  = i_len - step;
+		    for (j = 0; j < bound; j++) {
+			i = j + step + 1;
+			if (ta_popC[j] > ta_popC[i-1]) {
+			    for (k = 0; k < i_D; k++) 
+				tempP[k] = ta_popP.at(i-1, k);
+			    tempC = ta_popC[i-1];
+			    for (k = 0; k < i_D; k++) 
+				ta_popP.at(i-1,k) = ta_popP.at(j,k);
+			    ta_popC[i-1] = ta_popC[j];
+			    for (k = 0; k < i_D; k++) 
+				ta_popP.at(j,k) = tempP[k];
+			    ta_popC[j] = tempC;
+			    done = 0; 
+			    /* if a swap has been made we are not finished yet */
+			}  /* if */
+		    }  /* for */
+		} while (!done);   /* while */
+	    } /*while (step > 1) */
+	    /* now the best NP are in first NP places in gta_pop, use them */
+	    for (i = 0; i < i_NP; i++) {
+		for (j = 0; j < i_D; j++) 
+		    ta_newP.at(i,j) = ta_popP.at(i,j);
+		ta_newC[i] = ta_popC[i];
+	    }
+	} /*i_bs_flag*/
 
-	    if( trace > 0 ) {
-		if( (i_iter % 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]);
-		    Rprintf("\n");
+	/* have selected NP mutants move on to next generation */
+	for (i = 0; i < i_NP; i++) {
+	    for (j = 0; j < i_D; j++) 
+		ta_oldP.at(i,j) = ta_newP.at(i,j);
+	    ta_oldC[i] = ta_newC[i];
+	}
+	/* check if the best stayed the same, if necessary */
+	if(i_check_winner)  {
+	    same = 1;
+	    for (j = 0; j < i_D; j++)
+		if (t_bestitP[j] != t_bestP[j]) {
+		    same = 0;
 		}
+	    if(same && i_iter > 1)  {
+		i_xav++;
+		/* if re-evaluation of winner */
+		tmp_best = evaluate(l_nfeval, t_bestP.memptr(), par, fcall, rho);
+		
+		/* 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);
+		else
+		    t_bestC = tmp_best;
 	    }
-	} /* end loop through generations */
+	    else {
+		i_xav = 1;
+	    }
+	}
+	for (j = 0; j < i_D; j++) 
+	    t_bestitP[j] = t_bestP[j];
+	t_bestitC = t_bestC;
 
+	if( trace > 0 ) {
+	    if( (i_iter % 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]);
+		Rprintf("\n");
+	    }
+	}
+    } /* end loop through generations */
+    
     /* last population */
     k = 0;
     for (i = 0; i < i_NP; i++) {
@@ -506,7 +501,6 @@
 	    k++;
 	}
     }
-
     *gi_iter = i_iter;
 
     PutRNGstate();



More information about the Rcpp-commits mailing list