Using RcppArmadillo for a matrix population model

Standard

I’ve had a busy winter and spring with training for the Birkie and then recovering from the Birkie, hence the few posts. One of the things I’ve been doing is teaching myself the Rcpp package in R. This package lets people easily use C++ code within R. In this post, I demonstrate how I used Rcpp and specifically the RcppArmadillo package to create a population model.

Matrix models are popular in ecology. These models are series of difference equations (i.e., discrete time). I was interested in coding a simple example for a two life-stage species with the projection matrix R.

To code this up in R, I use the following code

popModelR A = matrix(c(0, 2, 0.25, 0.5), nrow = 2, byrow = TRUE),
P0 = c(0, 50)){
P = matrix(0, nrow = dim(A)[1], ncol = nYears + 1)
P[, 1] = P0
for( t in 1:nYears){
P[ , t + 1] = A  %*% P[, t]
}
return(P)
}

nYears = 10
A = matrix(c(0, 2, 0.25, 0.5),
nrow = 2, byrow = TRUE)
P0 = matrix(c(0, 50), nrow = 2)
P0
popModelROut <- popModelR(nYears = 10,
A = A,
P0 = P0)

Obviously, this simple model runs fast, but how would one code this with Rcpp? In order to get matrix operators, I needed to use the RcppArmadillo package, so my code looks like this:

library(“inline”)
library(“Rcpp”)
library(“RcppArmadillo”)

src1 <- ‘
int nYearsX = Rcpp::as<int>(nYears);
arma::mat P0X = Rcpp::as<arma::mat>(P0);
arma::mat AX  = Rcpp::as<arma::mat>(A);
arma::mat PX(AX.n_cols, nYearsX + 1);
PX.col(0) = P0X;

for(int t = 0; t < nYearsX; t++) {
PX.col(t + 1) =  AX * PX.col(t);
}

return Rcpp::wrap(PX);

popModelRcpp <- cxxfunction(signature(nYears = “integer”,
A = “matrix”,
P0 = “matrix”),
src1, plugin = “RcppArmadillo”)

popModelRcpp(nYears, A, P0)

Now, to compare the two functions, I use the benchmark package and run the model for 100 simulated years:

library(rbenchmark)
nYears = 100
res popModelR(nYears, A, P0),
columns = c(“test”, “replications”, “elapsed”,
“relative”, “user.self”, “sys.self”),
order = “relative”,
replications = 1000)
print(res)


test replications elapsed relative user.self sys.self
1 popModelRcpp(nYears, A, P0) 1000 0.02 1.00 0.00 0.00
2 popModelR(nYears, A, P0) 1000 0.53 29.61 0.64 0.00

The Rcpp code is almost 30 times quicker than the base code in R!

6 tips for a new LaTeX user

Standard
Recently a coworker started using LaTeX and asked for some tips. Here’s my 6 tips for starting to use LaTeX:
  1.  Start what you finish (i.e., close environments or else you get errors or weird bugs), for example \begin{equation} needs an \end{equation}
  2. Every document needs 3 things: \documentclass{<class>}, \begin{document}, and \end{document}
  3. For equations, use \( \) inline and \begin{eqnarray} \end{eqnarray} equations. \\ creates a new line. Use \\ \nonumber to continue on an equation and &=& to space multiple equations for example:
    <code>
    \begin{eqnarray}
    a &=& b +c \\
    a & =& b \\ \nonumber
     & & c
    \end{eqnarray}
    </code>
     Gives you something like:
    a = b +c   (1)
    a = b

          c        (2)

  4. Bib files are your friend for citations. Use Google Scholar to populate new citations.
  5.  \textit{My italics text}, \textbf{my bold text}, should get most of your formatting. Do NOT use the depreciated {\bf bold} or {\it italics} style. (cf http://tex.stackexchange.com/questions/41681/correct-way-to-bold-italicize-text for more details on the second point)
  6. {} can be very helpful, especially for complicated math functions, when order of operation is important. For example, \sigma_{\pi^2}^{2}

 

ArcPy tips and observations

Standard

When I started using ArcGIS a few weeks ago, I quickly became frustrated when operations failed after pointing and clicking several times. This motivated me start scripting in ArcPy. Now, I only need to load the script file once to make it crash. Thus, I am crashing more efficiently :). Overall, my experience with ArcPy has been alright. I discover two blog posts about hating ArcGIS 10 and ArcGIS sucking that are well worth the read. Despite the titles, the author actually likes ArcGIS. I understand why. ArcPy and ArcGIS are powerful. Plus, the documentation is usually well done and easy to find online (although, including data and reproducible examples like R’s help files would be a great addition).  However, ArcPy does have some problems/idiosyncrasies. Here’s what I have discovered:

1) Non-basic error message are often cryptic and some require hours of using Google. The archives on GIS Stack Exchange has become a new favorite page of mine. As an example, I got this error message (I removed and shortened my file paths in the code snip):

Layer (Folder) exists, please delete output if you wish to recreate layer
done with extraction for us_asp
0:00:36.046000
Could not create tile. name: z:\my documents\...\slopeclipped, adepth: 32, type: 1, iomode: 3, version: 3, compression: 0, eval case: 1 (GetOpenTileChannel)
Traceback (most recent call last):
File "", line 1, in
File "z:/My Documents/.../PythonFile.py", line 77, in
mF.RasterClipper(Slope, grid, SlopeClipped, arcpy.env.workspace)
File "myFunctions.py", line 37, in RasterClipper
arcpy.Clip_management(raster, clipExt, out, "", "ClippingGeometry")
File "C:\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 12015, in Clip
raise e
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to copy raster dataset
Failed to execute (Clip).

This error leads me to my second tip.

2) The workspace should not include a space. Sometimes, but not all the time, this causes problems. This took me a while to figure out and I only found this by chance online when looking up another problem. For example, this code is okay:
arcpy.env.workspace = "C:/Folder/"

is okay, but
arcpy.env.workspace = "C:/My Documents/"

caused my example problem.

3) Be prepared to use libraries such as NymPy. ArcPy does not always play well with USA wide data. This seems to be the case for ArcPy’s wrapper functions more so than the base functions (My guess is the best performing functions are GDAL functions). I plan on posting this function once I have written it and tested it.

4) “Python Scripting for ArcGIS” is still a great reference!

I plan to post more of my project as I work through it.

Learning Python

Standard

Over the past year, I have been teaching myself Python. I am writing this blog post so I remember how I wish I had learned Python in case I ever need to teach anybody else. My desire to learn Python has been driven by several factors: 1) my laziness and disdain for GUIs (i.e., I hate pointing and clicking); 2) my need to solve complex problems that GUIs would not solve (i.e., GUIs could not solve my problem); and 3) my curiosity to learn new things.

Python is a powerful computer language. I want to use Python for scripting and scientific computing, but many other uses of Python exist. A great way to learn Python and some computer science theory is with the book “Think Python: How to Think Like a Computer Scientist”. The book was originally written by the author for his college class using the Java language. However, he published the book with Creative Commons License and a high school teacher adapted the book for Python! Score one for open source freedom.

After working through this book, I went to the “Invent with Python” books, which are fun (but not necessarily that useful for my goal…). The first book, “Invent Your Own Computer Games with Python” covers basic game designs that are terminal based. The second book, “Making games with Python and Pygame” includes graphics. These books are good because they are simple and fun. If I was simply learning Python for fun, I would start with these books rather than “Think Python”. These books also gave me more experience coding (It’s been said the first 10,000 lines of coding are learning, so at least these gave me a fun couple of thousand), so perhaps the sugar coating of game development worked with me as well.

After these books, dive into subject area books. For me, these were GIS books. “Python Scripting for ArcGIS” does a great job introducing and covering ESRI’s flagship product, ArcGIS. The only downside is that the script is closed source. I decided to stop trying to read source code when a .py file including a warning about being a trade secret. “Learning Geospatial Analysis with Python” is an introduction to advanced (but open source) libraries such as Shapely and GDAL. However, the book is a bit pricey for what you get. But, the book is easier to read than the author’s blog that contains most of the same material. Also note that the ebook kills Pythons tab spacing (Python uses white space such as tabs to define functions and code). Also, I was initially overwhelmed because I did not know GIS or Python my first time reading it. My co-workers copy of “The Geospatial Desktop” provided a great refresher to my GIS skills!

Finally, Learning Python, 5th ed. is a tome, but very through. My own course of study began with the games books. Then I tried “Learning Geospatial Analysis with Python”. Along the way, I bought “Learning Python” and then discovered “Think Python”. I discovered the ArcGIS book only when I needed to script AcrPy for work. Not an optimal course of study, but needing to script ArcGIS has been a great crash course in Python for me. As an added bonus, Python has also made me a better R programmer! Horray for learning.