[Rcpp-devel] Problems wit rlnorm

Hi David,

| I?ve made a lot of progress working through some simple issues, but this one
| effects the cornerstone of my project.
| 
| Below is some sample code that shows my efforts.
| 
| I can set the seed programmatically with any integer variable in the place of
| the 20 in  ?seed?=20 using the environment method.

You can do that more simply at the R level. So I would just do

    set.seed(20)
    fun()         # calling the function you defined
 
| The language method could only make a fixed call.  I still can?t figure out
| what the RNGScope could possibly do.

Can you possibly go back and read the list archives? That was discussed at
some length when we did most of the sugar work -- in essence it just brings
the sane state of the R RNGs to our C++ projects, and then returns it. All
automatically witin the scope.
 
| I reviewed the discussion at:
| 
| http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2010-September/
| 001078.html
| 
| This made me realize that I should go ahead and code a simple putcol function
| loop when I need to populate a column in a Matrix ?array?.  (I had spent some
| time looking for such a defined method.)

Yes, you should also be able to assign a column directly. There are other examples.
 
| An array or Matrix class would be nice, but for now I seem to be able to find
| work-arounds.
| 
| Now I am calling the r-prefixed distributions with success, but this same code
| will fail on compile if the rlnorm() function is called.
| 
| src <- '               
| Rcpp::RNGScope Scope;  
| Environment base("package:base");      
| Function SetSeed = base["set.seed"];   
| int NumRands = 5;      
| int NumTrials = 3;     
| Rcpp::NumericVector RandCol(NumRands); 
|  Rcpp::NumericMatrix RandVals(NumRands, NumTrials*2);  
| SetSeed(Named("seed",20) );    
| RandCol=rnorm(NumRands,2,1);   
| int j=3;               
| for(int i=0; i<NumRands; i++) {
| RandVals(i,j) = RandCol(i);    
| }              
| return RandVals;       
| '              
|  fun <- cxxfunction(signature(),       
|  src, plugin = "Rcpp") 

I would write this as follows (forgetting for a moment that you do not need
to loop to assign the NumRands values into the matrix): 

src <- '
  Rcpp::RNGScope Scope;
  int NumRands = 5;
  int NumTrials = 3;
  Rcpp::NumericVector RandCol(NumRands);
  Rcpp::NumericMatrix RandVals(NumRands, NumTrials*2);
  RandCol = rnorm(NumRands,2,1);
  int j=3;
  for (int i=0; i<NumRands; i++) {
    RandVals(i,j) = RandCol(i);
  }
  return RandVals;
'
fun <- cxxfunction(signature(), body = src, plugin = "Rcpp", include="#include <Rcpp/stats/random/rlnorm.h>")
set.seed(20)
fun()

so no need for the Language or Enviornment objects.
                
| compile error if rlnorm is called instead of rnorm:    
| "file3a7d3a6f.cpp:38:28: error: 'rlnorm' was not declared in this scope"       
| 
| I looked at the source files and could not figure out why this does not work. 

I am a little stumped too. I also failed to get it to work. I will take
another look later --- this may in fact be a bug.
 
| Fortunately, the environment facility enabled me to have ?work around? success.
| 
| Adding these lines to the front of the  code solves the problem for now.
| 
| Environment stats("package:stats");    
| Function rlnorm = stats["rlnorm"];     

Yes, very good trick to get it from R -- but you should also be able to
access the rlnorm functions defined via include/Rcpp/stats/random/rlnorm.h

Dirk

Gauss once played himself in a zero-sum game and won $50.
                      -- #11 at http://www.gaussfacts.com

David,

| I am a little stumped too. I also failed to get it to work. I will take
| another look later --- this may in fact be a bug.

It _was_ a bug and seemingly nobody but you ever tried rnlorm :)  

There was a silly typo bug inasmuch as the header guard used _norm_ not
_lnorm_ and with the other term already defined, this header file was never
included.  We also need a few more fixed making exp() explicitly call the C
math library one, and correcting a typo.  

If you are comfortable applying a diff, try this one from rev3038.

Cheers,  Dirk


Index: inst/include/Rcpp/stats/random/rlnorm.h
===================================================================
--- inst/include/Rcpp/stats/random/rlnorm.h	(revision 3029)
+++ inst/include/Rcpp/stats/random/rlnorm.h	(working copy)
@@ -19,8 +19,8 @@
 // You should have received a copy of the GNU General Public License
 // along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef Rcpp__stats__random_norm_h
-#define Rcpp__stats__random_norm_h
+#ifndef Rcpp__stats__random_lnorm_h
+#define Rcpp__stats__random_lnorm_h
 
 namespace Rcpp {
 	namespace stats {
@@ -32,7 +32,7 @@
 				meanlog(meanlog_), sdlog(sdlog_) {}
 	
 			inline double operator()() const {
-				return exp( meanlog + sdlog * ::norm_rand() ) ;
+				return ::exp( meanlog + sdlog * ::norm_rand() ) ;
 			}
 	
 		private:
@@ -47,7 +47,7 @@
 				meanlog(meanlog_) {}
 	
 			inline double operator()() const {
-				return exp( meanlog + ::norm_rand() ) ;
+				return ::exp( meanlog + ::norm_rand() ) ;
 			}
 	
 		private:
@@ -57,10 +57,10 @@
 		class LNormGenerator_0 : public Generator<false,double> {
 		public:
 	
-			LNormGenerator_1( ) {}
+			LNormGenerator_0( ) {}
 	
 			inline double operator()() const {
-				return exp(::norm_rand() ) ;
+				return ::exp(::norm_rand() ) ;
 			}
 	
 		} ;
@@ -75,7 +75,7 @@
 			// TODO: R also throws a warning in that case, should we ?
 			return NumericVector( n, R_NaN ) ;
 		}  else if (sdlog == 0. || !R_FINITE(meanlog)){
-			return NumericVector( n, exp( meanlog ) ) ;
+			return NumericVector( n, ::exp( meanlog ) ) ;
 		} else {
 			return NumericVector( n, stats::LNormGenerator( meanlog, sdlog ) ); 
 		}

Gauss once played himself in a zero-sum game and won $50.
                      -- #11 at http://www.gaussfacts.com