tag:blogger.com,1999:blog-46954611923141121432024-03-13T06:16:41.615-07:00Life in CodeWork-products and stray thoughts from the Land of Entrapmentxianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.comBlogger109125tag:blogger.com,1999:blog-4695461192314112143.post-33027825469828235732016-11-14T03:21:00.000-07:002016-12-05T01:17:33.023-07:00Playing with OpenCLI spent last week reading up on modern C++ developments, including some great essays from <a href="https://herbsutter.com/">Herb Sutter</a>. I was particularly struck by his prescient series on Moore's Law, <a href="http://www.gotw.ca/publications/concurrency-ddj.htm">The Free Lunch Is Over</a> and <a href="https://herbsutter.com/welcome-to-the-jungle/">Welcome to the Jungle</a>. The latter essay portrays all possible computer architectures on a 2D plane of CPU versus memory architecture. The axes are a bit tricky, but the general idea is that a platform at the "origin" is predictable and easy to program for, whereas things get trickier as you move up and/or right.<br />
<br />
<p><a href="https://herbsutter.files.wordpress.com/2012/01/image3.png"><img style="display:inline;float:center;margin:0 0 0 10px;" title="image" src="https://herbsutter.files.wordpress.com/2012/01/image_thumb3.png?w=320&h=160" alt="image" width="320" height="160" align="center" /></a><br />
<br />
<p>This figure deserves describes everything from cellphones and game consoles to super computers and communications satellites. It also got me wondering how hard a simple "hello world" OpenCL program would be to get running on my Intel-only laptop. Can I do this in a day, or perhaps just an evening?<br />
<br />
Personally, I don't want to fuss with hardware right now - I just want to see how the GPU/C++ pieces fit together. Conveniently, my laptop contains a low-end 24-core embedded GPU on its Broadwell chip. Running Debian, I was able to easily <a href="https://github.com/helmingstay/test_code/blob/master/opencl/README.md">download the requisite packages</a> and get started. I quickly discovered that consumer Intel GPUs, at present, do *not* <a href="http://arrayfire.com/opencl-on-intel-hd-iris-graphics-on-linux/">support <code>double</code> anything</a>, making this a less-than-ideal test-bed for scientific programming.<br />
<br />
My test case was inspired by a performance issue that I ran into in my work. In a scientific simulation program that I use & help develop, Valgrind revealed that <code>pow(double, double)</code> was taking fully half of the total computational time. Poking around a bit, I see that <code>pow()</code> and <code>log</code> really are <a href="https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Elementary_functions">quite complex to compute</a>, particularly for doubles (since the total effort is a function of precision). With this in mind, I set up a simple example using both OpenCl and straight C++, and compared timings. Note - I strongly recommend using a sample size of greather than one to draw any conclusions with real-life consequences!<br />
<br />
In this example, the vanilla C++ is clean and easy to read, but is ~20x slower than the OpenCL version. Worried about the possibility of "unintended optimizations", I tried using a different kernel function. I used <code>float</code> for both examples to keep the total computational complexity the same. The speed results remained, but the new test revealed different answers. To the best of my understanding, this highlights differences in precision-sensitive operations between the OpenCL and stdc++ platforms. This is a pretty tricky area - just know that it's something to keep an eye on if you require perfect concordance between platforms.<br />
<br />
EDIT: I also added an example using <a href="http://boostorg.github.io/compute">Boost.Compute</a> today, which brings the best of both the C++ and OpenCL worlds. <code>Boost.Compute</code> has straighforward docs, and includes a nice <a href="http://boostorg.github.io/compute/BOOST_COMPUTE_CLOSURE.html">closure macro</a> that allows the direct incorporation of C++ code in kernel functions. The resulting code is *way* less verbose than vanilla OpenCL. The only downside is the extra requirement. That and some *very* noisy compiler warnings.<br />
<br />
Here's a full example that <a href="https://github.com/helmingstay/test_code/tree/master/opencl">can be found in my github test code repo</a> (boost not shown, timings comparable with OpenCL):<br />
<pre>$make; time ./opencl; time ./straight
g++ -Wall -std=c++11 -lOpenCL -o opencl opencl.cpp
g++ -Wall -std=c++11 -o straight straight.cpp
Using platform: Intel Gen OCL Driver
Using device: Intel(R) HD Graphics 5500 BroadWell U-Processor GT2
result:
200 201 202 203 204 205 206 207 208 209
99990 99991 99992 99993 99994 99995 99996 99997 99998 99999
./opencl 0.20s user 0.09s system 97% cpu 0.295 total
result:
200 201 202 203 204 205 206 207 208 209
99990 99991 99992 99993 99994 99995 99996 99997 99998 99999
./straight 6.03s user 0.00s system 99% cpu 6.032 total
</pre><br />
Hopefully this example helps you get started experimenting with GPU computing. As <a href="https://herbsutter.com/">Herb Sutter</a> points out, we can expect more and greater hardware parallelism in the near future. Discrete GPUs are now commonly used in scientific computing, and Intel is now selling a massively-multicore add-on card, <a href="https://en.wikipedia.org/wiki/Xeon_Phi">the Xeon Phi processor</a>. Finally, Floating point precision <a href="http://arrayfire.com/explaining-fp64-performance-on-gpus/">remains an interesting question</a> to keep an eye on in this domain.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-85251780680105143412016-02-15T22:06:00.000-07:002016-02-17T23:08:00.361-07:00Shiny on Webfaction: VPS installion without root<p>I've been using <a href="https://www.webfaction.com/?aid=88752">Webfaction (plug)</a> as an inexpensive managed VPN. Part of me wants VPS root access, but I'm mostly happy to leave the administrative details to others. Webfaction seems to be a good example of a common VPS plan: user-only access in a rich development environment. Compilers, <code>zsh</code>, and even <code>tmux</code> are available from the shell, making this a very comfortable dev environment overall.
<p>Most times root doesn't matter, but sometimes it complicates new software installs. I've been looking forwards to testing R's webapp package <a href="http://shiny.rstudio.com/">Shiny</a>, but all of the docs assume root access (and some even state that it's required). I set off without knowing if this would work, attempting to see how far I could get. What follows is a (hopefully) reproducible account of a user-land install of R & Shiny via ssh on a Webfaction slice. To the best of my knowledge, this requires only standard development tools, and so should(??) work.
<p>In the following I use [tab] to indicate hitting tab key for auto-completion. The VPS login username is [user]. [edit] means call your editor of choice (vim, emacs, or, god forbid, nano). This assumes you are using bash (which seems to be the default shell on most VPNs).
<h4> Prepare the build environment</h4>
<pre>
## ssh to webhost
## make directories, set paths, etc
## source build dir
mkdir ~/src
## software install dir
mkdir ~/local
## personal content dir
CONTENTDIR=~/var
mkdir $CONTENTDIR
## some hosts have /tmp set noexec?
mkdir src/tmp
## Install software here
INSTPREFIX=$HOME/local
## set paths:
##
echo 'export PATH=$PATH:~/local/bin:~/local/shiny-server/bin' >> ~/.bashrc
echo 'export TMPDIR=$HOME/src/tmp' >>~/.bashrc
## check that all is well
[edit] ~/.bashrc
## update env
. .bashrc
</pre>
[Ref: <a href="http://mazamascience.com/WorkingWithData/?p=1185">temp dir and R packages</a>]
<h4>Install R from source: fast and (mostly) easy</h4>
<pre>
cd ~/src
wget http://cran.us.r-project.org/src/base/R-3/R-3.2.3.tar.gz
tar xzf R-3.2.3.tar.gz
cd R-[tab]
./configure --prefix=$INSTPREFIX
## missing library, search and add directory
CPPFLAGS=/usr/lib/jvm/java/include/ make
make install
cd ~
</pre>
<h4>Prep R environment</h4>
<pre>
## The following commands are in R
install.packages(c('shiny', 'rmarkdown'))
</pre>
<pre>
## From the shell:
## on a headless / no-X11 box, need cairo for png
echo "options(bitmapType='cairo')" >> ~/.Rprofile
## check that all is well
[edit] ~/.Rprofile
</pre>
[Ref: <a href="http://stackoverflow.com/questions/13067751/how-to-run-r-scripts-on-servers-without-x11">R png without X11</a>]
<h4>Install cmake (if needed)</h4>
<pre>
## first install cmake - skip if's already available
`which cmake`
## nothing? continue
## NOTE - I'm using the source tarball here, not binaries
wget https://cmake.org/files/v3.4/cmake-3.4.3.tar.gz
tar xzf cmake-[tab]
cd cmake-[tab]
./configure --prefix=$INSTPREFIX
gmake
make install
</pre>
<h4>Install Shiny Server</h4>
<pre>
## From shell
cd ~/src
git clone https://github.com/rstudio/shiny-server.git
cd shiny-server
cmake -DCMAKE_INSTALL_PREFIX=$INSTPREFIX
make
## "make install" Complains about no build dir
## I'm not sure what happens here, but this seems to work
PYTHON=`which python`
mkdir build
./bin/npm --python="$PYTHON" rebuild
./bin/node ./ext/node/lib/node_modules/npm/node_modules/node-gyp/bin/node-gyp.js --python="$PYTHON" rebuild
make install
</pre>
[Ref: <a href="https://github.com/rstudio/shiny-server/wiki/Building-Shiny-Server-from-Source">shiny build docs</a>]
<h4>Configure Shiny Server</h4>
All of the <a href="https://rstudio.github.io/shiny-server/latest/#server-management">Shiny Server docs</a> assume the config file is located in <code>/etc/</code>, which I don't have access to. There's _zero_ documentation on running shiny, nor does running <code>shiny-server -h</code> or <code>shiny-server --help</code> provide any indication. Trial and error and reading source code on github finally leads to <code>shiny-server path-to-config-file</code>. So, let's make a shiny site!
<pre>
## Nest content in ~/var
mkdir $CONTENTDIR/shiny
cp -rp ~/src/shiny-server/samples $CONTENTDIR/shiny/apps
mkdir $CONTENTDIR/shiny/logs
## copy the packaged settings template to the content dir
cp ~/src/shiny-server/config/default.config $CONTENTDIR/shiny/server.conf
[edit] $CONTENTDIR/shiny/server.conf
##
## server.conf content follows:
run_as [user];
## leave location as-is
## substitute var with $CONTENTDIR if needed
site_dir /home/[user]/var/shiny/apps;
log_dir /home/[user]/var/shiny/logs;
## save file
## back at shell, run shiny, put in background
shiny-server ~/var/shiny/server.conf &
</pre>
[Ref: <a href="http://rstudio.github.io/shiny-server/latest/#r-markdown">Shiny-server docs</a>]
<h4>Testing</h4>
Shiny should give messages about <code>Starting listener on 0.0.0.0:3838</code>. First up, let's use ssh to connect remote port 3838 to a local port. This allows local testing without deployment. As an aside, if you're not using <code>~/.ssh/config</code> on a local machine to manage keys and hostname shortcuts, you should!
<pre>
## on local machine:
ssh -nNT -L 9000:127.0.0.1:3838 [user]@webhost
</pre>
Now, if all went well, you should be able to navigate to the welcome page via browser on local machine:
<br><code>http://127.0.0.1:9000</code>
<p>Once shiny is working, don't forget to take a look at your logs:
<br><code>ls -alh $CONTENTDIR/shiny/logs</code>
<p>I had trouble with the packaged <code>rmd</code> example app (which renders a <code>.Rmd</code> file). Reading logs showed install issues with pandoc, and I had to manually fiddle with the links:
<pre>ln -s $INSTPREFIX/shiny-server/ext/pandoc/static/pandoc $INSTPREFIX/shiny-server/ext/pandoc/</pre>
[Ref: <a href="http://blog.trackets.com/2014/05/17/ssh-tunnel-local-and-remote-port-forwarding-explained-with-examples.html">port forwarding</a>]
<h4>Wrap-up</h4>
For a full production environment, you would want a process monitor to keep <code>shiny-server</code> running, as well a public-facing server. See your webhost's documentation for process monitors. More details on <a href="https://support.rstudio.com/hc/en-us/articles/213733868-Running-Shiny-Server-with-a-Proxy">shiny-server and apache are here</a> (I haven't tried these proxy methods).
<p>Finally, a more conventional approach using root access on a VPS (such as DigitalOcean) is <a href="http://deanattali.com/2015/05/09/setup-rstudio-shiny-server-digital-ocean/">available here</a>.
<h4>Update - 17 Feb 2016: Deployment Logistics</h4>
After a day of kicking the tires, I'm happy to report Shiny-server is working well on Webfaction in production mode. Two points:
<p><b>Making a webapp.</b> In the Webfaction control panel, I added a custom application. In the following, substitute [appname] for the value entered in the Name field. For <code>App category</code> I selected "Websockets", and then clicked "Save". Copy the port number. Edit the <code>server.conf</code> file from above, replacing the number in <code>listen 3838;</code> with the port number copied from Webfaction. Finally, create a website, add a name (can be the same as [appname] from above), and a domain. It typically takes a few minutes for DNS changes to propagate.
<p>The above steps creates a directory named <code>$HOME/webapps/[appname]</code>. I placed the <code>server.conf</code> file here, created app and log directories, and then updated server.conf to reflect the new locations:
<pre>
## Create the following directories
## add these paths to server.conf,
## and don't forget the trailing ;
mkdir $HOME/[appname]/logs
## shiny app files go here:
mkdir $HOME/[appname]/app
</pre>
<p>[Ref: <a href="https://docs.webfaction.com/software/custom.html">Webfaction custom application</a>]
<br>[Ref: <a href="https://docs.webfaction.com/user-guide/websites.html#websites-create">Webfaction Applications and Websites</a>]
<p><b>Running the server.</b> Shiny-server will use a PID file, which makes job-spawning a simple shell script + cron job. If shiny-server is already running, it will recognize the PID file and not start another process. I made the following script:
<pre>
#!/bin/sh
## executable shell script names $HOME/bin/my.shiny.sh
## make sure to run: chmod +x $HOME/bin/my.shiny.sh
APPROOT=$HOME/webapps/[appname]
PIDFN=$APPROOT/shiny-server.pid
## using full path
$HOME/local/shiny-server/bin/shiny-server $APPROOT/server.conf --pidfile=$PIDFN>> $APPROOT/logs/server.log 2>&1 &
</pre>
Now run <code>crontab -e</code> and add an entry for the script (above):
<pre>
## try once an hour, on the 10th minute of the hour
10 * * * * /home/[user]/bin/my.shiny.sh
</pre>
Finally, take a look at memory usage. If you exceed memory limits, Webfaction automatically kills everything. And R's memory use grows with more connections (which themselves persist, because <a href="https://en.wikipedia.org/wiki/WebSocket">websockets</a>). Webfaction distributes a <a href="https://wiki.webfaction.com/attachment/wiki/MiscellaneousFiles/memory_usage.py">nice python script</a> that shows per-process and total memory usage.
<p>[Ref: <a href="https://github.com/rstudio/shiny-server/blob/master/config/systemd/shiny-server.service">shiny-server systemd script</a> (shows commandline usage)]
<br>[Ref: <a href="https://docs.webfaction.com/software/general.html#scheduling-tasks-with-cron">Webfaction cron</a>]
<p>I should point out that I like <a href="https://www.webfaction.com/?aid=88752">Webfaction (plug)</a> well enough to pay them money. Their intro plan is $10/month for 1GB RAM + 100GB full SSD, with a 1-month free trial. I like that the webfaction user-base is big enough that lots of my questions are already answered, but small enough that staff actually answer new questions.
<p>I've done my best to document exactly what I did, but I'm sure there are typos. Let me know if you encounter any issues!
xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com8tag:blogger.com,1999:blog-4695461192314112143.post-88864895886398060322014-11-14T01:53:00.000-07:002014-11-14T02:11:10.390-07:00Numerical Simulations and Data passing: C++, Python, and Protocol Buffers <h2>Problem statement & Requirements</h2>
I'm working with a complex C++ simulation that requires a large number of user-specified parameters. Both speed and readability are important. I'd like to define all possible parameters in one (and only one) place, and include sensible defaults that can be easily over-ridden. Finally, intelligent type-handling would be nice.
For convenience, I decided to wrap the C++ simulation in python setup/glue code. Python is a logical choice here as the "available everywhere" glue language that has nice standard libraries.
<h2>Available libraries</h2>
There aren't many data-passing options that work with both C++ and python. Libconfig, JSON, XML, and Google Protocol Buffers (PB) appear to be the only reasonable options. Here's my thoughts on the first three:
<ul>
<li><a href="http://www.hyperrealm.com/libconfig/">Libconfig</a>: Nice clean library, good language support. The big downside is that data structures must be defined both in a data file and in code - e.g. data is "moved" from a file into C++ variables. I feel like libconfig is best for a small number of complex variables, like lists and vectors.
<li><a href=""http://www.thomaswhitton.com/blog/2013/06/28/json-c-plus-plus-examples/>JSON</a>: no clear standard C++ library, library docs so-so, speed complaints from some?
<li>XML: Massive overkill.
</ul>
That leaves PB, which has <a href="https://developers.google.com/protocol-buffers/docs/proto">nice docs</a> for both C++ and python. All the variables, along with their types and defaults, are defined in a <tt>.proto</tt> file. The <tt>protoc</tt> tool auto-generates python and C++ code from the .proto file. By adding it to my Makefile, C++ classes are autogenerated at compile time. This makes for fast and readable C++ code - like using a named dict, but without the speed costs.
<h2> Solution / Workflow </h2>
I'm using python to read user-supplied values into a set of PB messages, and then serializing the messages to files. C++ then reads the messages from those files at runtime. A python script run by <tt>make</tt> synchronizes the locations of files between python and C++.
I also want to process commandline options for my python wrapper script. Happily, I can hand a PB message to <a href="https://docs.python.org/3.1/library/optparse.html#optparse-reference-guide">python's parser.parse_args()</a> and have it <a href="http://comments.gmane.org/gmane.comp.lib.protocol-buffers.general/2040
">set PB message attributes with setattr()</a>.
The last python step (aside from writing the message to disk) is reading "variable,value" pairs from a .csv file. If a variable has already been set by parse_args, I skip it: the commandline values override .csv file values.
<h2>Summary</h2>
Overall, PB makes a very nice data coupler between an interpreted language like python and a compiled language like C++. Python excels at text processing and is easy to prototype, while C++ is fast and beautiful.
PB has a few side-benefits. On the C++ side, it provides some natural namespace encapsulation to manage variable-explosion. Runtime inspection with gdb is easy enough.
Finally, storing all the options values used to run each simulation in a standard-format file is handy - it allows tests to re-run the simulation with exactly the same inputs.
<h2>Python Snippets</h2>
<pre>
def main():
## initialize protobuf, fill with ParseArgs
setupSim = ProtoBufInput_pb2.setupSim()
setupSim = ParseArgs(sys.argv[1:], setupSim)
prepInput(setupSim)
RunSim()
def ParseArgs(argv, setupSim):
parser = OptionParser(usage="wrapper.py [options]\nNote: commandline args over-ride values in files.", version=setupSim.version)
## these must be valid protocol buffer fields
parser.add_option("-t", "--test", dest="testCLI",
action='store_true', help="Run test suite")
parser.add_option("-d", "--days", metavar='N',
dest="number_of_days",
type='int', help="Number of days to simulate")
## parse!
(setupSim, args) = parser.parse_args(argv, values=setupSim)
return(setupSim)
def prepInput(setupSim):
## options from ParseArgs
inhandle = open(setupSim.file_options, 'r')
outhandle = open(ProtoDataFiles.PbFile_setupSim, 'wb')
reader = csv.reader(inhandle, delimiter=',')
header = reader.next()
if not (header == ['variable','value']):
raise Exception('Incorrect header format')
for row in reader:
## skip comments, check for 2 fields per row
if (row[0][0] == '#'):
continue
if not (len(row) == 2):
raise Exception('Problem with value pair: %s' % row)
## pack the message using text representation
msgText = '%s : %s' % (row[0], row[1])
if setupSim.HasField(row[0]):
print("Skipping config file, keeping commandline value: %s, %s" % (row[0], getattr(setupSim,row[0])))
continue
setupSim = Merge(msgText, setupSim)
## write out to file for C++ to read
outhandle.write(setupSim.SerializeToString())
outhandle.close()
def RunSim():
subprocess.Popen("./sim").communicate()
if __name__ == "__main__":
main()
</pre>
<h2>C++ Code Snippets</h2>
<pre>
//PbRead.h
#include <iostream>
#include <fstream>
#include <string>
#include <stdexcept>
#include "proto/ProtoBufInput.pb.h"
template <typename Type>
void PbRead(Type &msg, const char *filename){
std::fstream infile(filename, std::ios::in | std::ios::binary);
if (!infile) {
throw std::runtime_error("Setup message file not found");
} else if (!msg.ParseFromIstream(&infile)) {
throw std::runtime_error("Parse error in message file");
}
}
// sim.cpp
#include "PbRead.h"
#include "ProtoDataFiles.h"
// protocol buffers get passed around, are globals
ProtoBufInput::setupSim PbSetupSim;
int main(int argc,char **argv)
{
GOOGLE_PROTOBUF_VERIFY_VERSION;
PbSetupSim.set_init(true);
// #define PbFile_setupSim "filename" in ProtoDataFiles.h, written by make
PbRead(PbSetupSim, PbFile_setupSim);
//...
if (PbSetupSim.test_2()){
//...
}
}
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-20003429217128502972014-07-03T19:38:00.000-07:002014-07-03T20:14:46.323-07:00Efficient Ragged Arrays in R and Rcpp<h2>When is R Slow, and Why?</h2>
Computational speed is a common complaint lodged against R.
Some recent posts on <a href="http://r-bloggers.com">r-bloggers.com</a> have compared
the speed of R with some other programming languages [<a href="#ref1">1</a>], and showed the
favorable impact of the new compiler package on run-times [<a href="#ref2">2</a>]. I and
others have written about using Rcpp to easily write C++ functions to speed-up bottlenecks
in R [<a href="#ref3">3</a>,<a href="#ref4">4</a>]. With the new Rcpp attributes framework,
writing fully vectorized C++ functions and incorporating them in R code is now very easy [<a
href="#ref5">5</a>].
<br>
<br>
On a day-to-day basis, though, R's performance is largely a function of coding style.
R allows novices users to write horribly inefficient code [<a href="#ref6">6</a>]
that produces the correct answer (eventually). Yet by failing to utilize vectorization
and pre-allocation of data structures, naive R code can be many orders of magnitude slower than
need be. <a href="https://stat.ethz.ch/mailman/listinfo/r-help">R-help</a> is littered with the
tears of novices, and there's even a (fantastic)
parody of Dante's Inferno outlining the common "Deadly Sins of R" [<a href="#ref7">7</a>].
<br>
<h2>Problem Statement: Appending to Ragged Arrays</h2>
I recently stumbled onto an interesting code optimization problem that I *didn't* have a quick solution for,
and that I'm sure others
have encountered. What is the "R way" to vectorize computations on ragged array?
One example of a ragged array is
a list of vectors that have varying and different lengths. Say you need to dynamically grow
many vectors by varying lengths
over the course of a stochastic simulation. Using a simple tool like <code>lapply</code>, the entire data
structure will be allocated anew with every assignment. This problem is briefly touched
on in the official
<a href="http://cran.r-project.org/doc/manuals/R-intro.html#The-function-tapply_0028_0029-and-ragged-arrays">Introduction
to R</a> documentation, which simply notes that "when the subclass sizes [e.g. vector sizes] are all the same the indexing
may be done implicitly and much more efficiently". But what if you're data *isn't* rectangular? How
might one intelligently vectorize a ragged array to prevent (sloooow) memory allocations at every step?
<br>
<br>
The obvious answer is to pre-allocate a rectangular matrix (or array) that is larger than the maximum
possible vector length, and store each vector as a row (or column?) in the matrix. Now we can use matrix
assignment, and for each vector track the index of the start of free space.
If we try to write past the end of the matrix, R emits the appropriate error.
This method requires some book-keeping on our part. One nice addition would be
an S4 class with slots for the data matrix and the vector of free-space
indices, as well as methods to dynamically expand the data matrix and validate the object.
As an aside, this solution is essentially the inverse of a sparse matrix. Sparse matrices use
much less memory at the expense of slower access times [<a href="#ref8">8</a>].
Here, we're using more memory than is strictly needed to achieve much faster access times.
<br>
<br>
Is pre-allocation and book-keeping worth the trouble? <code>object.size(matrix(1.5, nrow=1e3,
ncol=1e3))</code> shows that a data structure of 1,000 vectors, each of length approximately 1,000,
occupies about 8Mb of memory. Let's say I resize this structure 1,000 times. Now I'm looking
at almost a gigabyte of memory allocations. Perhaps you're getting a sense of what a terrible idea
it is to *not* pre-allocate a frequently-resized ragged list?
<br>
<h2>Three Solutions and Some Tests</h2>
Using the above logic, I prototyped a solution as an R function, and then transcribed the result into
a C++ function (boundary checks are important in C++).
The result is three methods: a "naive" list-append method, an R method that uses matrix assignment,
and a final C++ method that modifies the pre-allocated matrix in-place.
In C++/Rcpp, functions can use pass-by-reference semantics [<a href="#ref9">9</a>],
which can have major speed advantages by
allowing functions to modify their arguments in-place. Full disclosure: pass-by-reference semantics requires
some caution on the user's part. Pass-by-reference is very different from R's "function
programming" semantics (pass-by-value, copy-on-modify), where side-effects are minimized and an explicit
assignment call is required to modify an object [<a href="#ref10">10</a>].
<br>
<br>
I added a unit test to ensure identical results between all three methods, and then used the fantastic <a
href="http://cran.r-project.org/web/packages/rbenchmark/index.html">rbenchmark</a> package to time each solution.
As expected, the naive method is laughably slow. By comparison, and perhaps counter-intuitively, the R and C++
pre-allocation methods are close in performance. Only with more iterations and larger
data structures does the C++ method really start to pull ahead. And by that time, the naive R
method takes *forever*.
<br>
<br>
Refactoring existing code to use the pre-allocated compound data structure (matrix plus indices)
is a more challenging exercise that's "left to the reader", as mathematics textbooks oft say.
<code>lapply()</code> is
conceptually simple, and is often fast *enough*. Some work is required to transcribe code from this
simpler style to use the "anti-sparse" matrix (and indices). There's a temptation to prototype a solution
using <code>lapply()</code>
and then "fix" things later. But if you're using ragged arrays and doing any heavy lifting
(large data structures, many iterations), the timings show that pre-allocation is more than worth the effort.
<br>
<h2>Code</h2>
Note: you can find the full code <a href="http://www.x14n.org/code/raggedArrayTest.R">here</a> and
<a href="http://www.x14n.org/code/helpers.cpp">here</a>.
<br><br>
Setup: two helper functions are used to generate ragged arrays via random
draws. First, draws from the negative binomial distribution determine the length of the each new vector (with a
minimum length of 1, <code>gen.lengths()</code>), and draws from the normal distribution fill each
vector with data (<code>gen.dat()</code>).
<pre>
## helper functions
gen.lengths <- function(particles, ntrials=3, prob=0.5) {
## a vector of random draws
pmax(1, rnbinom(particles, ntrials, prob))
}
gen.dat <- function(nelem, rate=1) {
## a list of vectors, vector i has length nelem[i]
## each vector is filled with random draws
lapply(nelem, rexp, rate=rate)
}
</pre>
<br>
Three solutions: a naive <code>lapply()</code> method, followed by pre-allocation
in R.
<pre>
## naive method
appendL <- function(new.dat, new.lengths, dat, dat.lengths) {
## grow dat by appending to list element i
## memory will be reallocated at each call
dat <- mapply( append, dat, new.dat )
## update lengths
dat.lengths <- dat.lengths + new.lengths
return(list(dat=dat, len=dat.lengths))
}
## dynamically append to preallocated matrix
## maintain a vector of the number of "filled" elements in each row
## emit error if overfilled
## R solution
appendR <- function(new.dat, new.lengths, dat, dat.lengths) {
## grow pre-allocated dat by inserting data in the correct place
for (irow in 1:length(new.dat)) {
## insert one vector at a time
## col indices for where to insert new.dat
cols.ii <- (dat.lengths[irow]+1):(dat.lengths[irow]+new.lengths[irow])
dat[irow, cols.ii] = new.dat[[irow]]
}
## update lengths
dat.lengths <- dat.lengths + new.lengths
return(list(dat=dat, len=dat.lengths))
}
</pre>
<br>
<br>
Next, the solution as a C++ function. This goes in a separate file that I'll call <code>helpers.cpp</code> (compiled below).
<pre>
#include <Rcpp.h>
using namespace Rcpp ;
// [[Rcpp::export]]
void appendRcpp( List fillVecs, NumericVector newLengths, NumericMatrix retmat, NumericVector retmatLengths) {
// "append" fill oldmat w/
// we will loop through rows, filling retmat in with the vectors in list
// then update retmat_size to index the next free
// newLenths isn't used, added for compatibility
NumericVector fillTmp;
int sizeOld, sizeAdd, sizeNew;
// pull out dimensions of matrix to fill
int nrow = retmat.nrow();
int ncol = retmat.ncol();
// check that dimensions match
if ( nrow != retmatLengths.size() || nrow != fillVecs.size()) {
throw std::range_error("In appendC(): dimension mismatch");
}
for (int ii = 0; ii<nrow; ii++){
// for each vector / row in retmat
// vector to add to retmat
fillTmp = fillVecs[ii];
// compute lengths
sizeOld = retmatLengths[ii];
sizeAdd = fillTmp.size();
sizeNew = sizeOld + sizeAdd;
// error checking - stop on overfill
if ( sizeNew >= ncol) {
throw std::range_error("In appendC(): exceeded max cols");
}
// iterator for row to fill
NumericMatrix::Row retRow = retmat(ii, _);
// fill row of return matrix, starting at first non-zero elem
std::copy( fillTmp.begin(), fillTmp.end(), retRow.begin() + sizeOld);
// update size of retmat
retmatLengths[ii] = sizeNew;
}
}
</pre>
<br>
<br>
Putting the pieces together: a unit test ensures the results of all three methods are identical, and a function that runs each
solution with identical data will be used for timing.
<pre>
## unit test
test.correct.append <- function(nrep, particles=1e3, max.cols=1e3, do.c=F) {
## list of empty vectors, fill with append
dat.list <- lapply(1:particles, function(x) numeric())
## preallocated matrix, fill rows from left to right
dat.r <- dat.c <- matrix(numeric(), nrow=particles, ncol=max.cols)
## length of each element/row
N.list <- N.r <- N.c <- rep(0, particles)
## repeat process, "appending" as we go
for (ii in 1:nrep) {
N.new <- gen.lengths(particles)
dat.new <- gen.dat(N.new)
## in R, list of vectors
tmp <- appendL(dat.new, N.new, dat.list, N.list)
## unpack, update
dat.list <- tmp$dat
N.list <- tmp$len
## in R, preallocate
tmp <- appendR(dat.new, N.new, dat.r, N.r)
## unpack, update
dat.r <- tmp$dat
N.r <- tmp$len
## as above for C, modify dat.c and N.c in place
appendRcpp(dat.new, N.new, dat.c, N.c)
}
## pull pre-allocated data back into list
dat.r.list <- apply(dat.r, 1, function(x) { x <- na.omit(x); attributes(x) <- NULL; x } )
## check that everything is
identical(dat.r, dat.c) && identical(N.r, N.c) &&
identical(dat.list, dat.r.list) && identical(N.list, N.r)
}
## timing function, test each method
test.time.append <- function(nrep, particles=1e2, max.cols=1e3, append.fun, do.update=T, seed=2) {
## object to modify
N.test <- rep(0, particles)
dat.test <- matrix(numeric(), nrow=particles, ncol=max.cols)
## speed is affected by size,
## so ensure that each run the same elements
set.seed(seed)
for (irep in 1:nrep) {
## generate draws
N.new <- gen.lengths(particles)
dat.new <- gen.dat(N.new)
## bind in using given method
tmp <- append.fun(dat.new, N.new, dat.test, N.test)
if(do.update) {
## skip update for C
dat.test <- tmp$dat
N.test <- tmp$len
}
}
}
</pre>
<br>
<br>
Finally, we run it all:
<pre>
library(rbenchmark)
## Obviously, Rcpp requires a C++ compiler
library(Rcpp)
## compilation, linking, and loading of the C++ function into R is done behind the scenes
sourceCpp("helpers.cpp")
## run unit test, verify both functions return identical results
is.identical <- test.correct.append(1e1, max.cols=1e3)
print(is.identical)
## test timings of each solution
test.nreps <- 10
test.ncols <- 1e3
timings = benchmark(
r=test.time.append(test.nreps, max.cols=test.ncols, append.fun=appendR),
c=test.time.append(test.nreps, max.cols=test.ncols, do.update=F, append.fun=appendRcpp),
list=test.time.append(test.nreps, max.cols=test.ncols, append.fun=appendL),
replications=10
)
## Just compare the two faster methods with larger data structures.
test.nreps <- 1e2
test.ncols <- 1e4
timings.fast = benchmark(
r=test.time.append(test.nreps, max.cols=test.ncols, append.fun=appendR),
c=test.time.append(test.nreps, max.cols=test.ncols, do.update=F, append.fun=appendRcpp),
replications=1e1
)
</pre>
<br>
<br>
Benchmark results show that the list-append method is 500 times slower than the improved
R method, and 1,000 times slower than the C++ method (<code>timings</code>).
As we move to larger data structures
(<code>timings.fast</code>), the advantage of modifying in-place with C++ rather than having to
explicitly assign the results quickly add up.
<pre>
> timings
test replications elapsed relative user.self sys.self
2 c 10 0.057 1.000 0.056 0.000
3 list 10 52.792 926.175 52.674 0.036
1 r 10 0.128 2.246 0.123 0.003
> timings.fast
test replications elapsed relative user.self sys.self
2 c 10 0.684 1.000 0.683 0.000
1 r 10 24.962 36.494 24.934 0.027
</pre>
<h2>References</h2>
<br><a name="ref1">[1]</a> <a href="http://www.r-bloggers.com/how-slow-is-r-really/">How slow is R really?</a>
<br><a name="ref2">[2]</a> <a href="http://www.r-bloggers.com/speeding-up-r-computations-pt-ii-compiling/">Speeding
up R computations Pt II: compiling</a>
<br><a name="ref3">[3]</a> <a
href="http://helmingstay.blogspot.com/2011/06/efficient-loops-in-r-complexity-versus.html">Efficient loops in R - the complexity versus speed trade-off</a>
<br><a name="ref4">[4]</a> <a href="http://dirk.eddelbuettel.com/blog/2011/09/08/">Faster (recursive) function calls: Another quick Rcpp case study</a>
<br><a name="ref5">[5]</a> <a href="http://dirk.eddelbuettel.com/blog/2012/11/20/">Rcpp attributes: A simple example 'making pi'</a>
<br><a name="ref6">[6]</a> <a href="http://stackoverflow.com/questions/2908822/speed-up-the-loop-operation-in-r">StackOverflow: Speed up the loop operation in R</a>
<br><a name="ref7">[7]</a> <a href="http://www.burns-stat.com/documents/books/the-r-inferno/">The R Inferno</a>
<br><a name="ref8">[8]</a> <a href="http://blog.custora.com/2013/05/sparse-matrix-formats-pros-and-cons/">Sparse matrix formats: pros and cons</a>
<br><a name="ref9">[9]</a> <a href="http://adv-r.had.co.nz/OO-essentials.html#rc">Advanced R: OO field guide: Reference Classes</a>
<br><a name="ref10">[10]</a> <a href="http://adv-r.had.co.nz/Functions.html#return-values">Advanced R: Functions: Return Values</a>
xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-41579363443935006882014-05-30T18:19:00.000-07:002014-05-30T18:24:33.641-07:00 Tools for Online Teaching<p>
Last semester (Fall 2014), I organized and taught an interdisciplinary,
collaborative class titled <a
href="http://probforsci.blogspot.com/">Probability for Scientists</a>.
Getting 4 separate teachers on the same page was a challenge, but as
scientists we're used to communicating over email, and CC'ing everyone
worked well enough. Throw in a shared dropbox folder, a shared google
calendar, and a weekly planning meeting, and instructor
collaboration went pretty smoothly.
<p>
It was more challenging to organize the class so that we could easily provide students with
up-to-date course information and supplemental material.
We ended up using blogger, which has some key benefits and
disadvantages. It was *really* easy to set up and add permissions for
multiple individuals. This allowed any of us to put up images (for
example, photos of the whiteboard at the end of class) and post
instructions. One downside we heard from students was an apparent
lack of organization. I attempted to organize the blog with
intelligent post labels, along with the "Labels widget" (which shows
all possible labels and the number of posts per label on the right
hand sidebar). I also included the "Posts by Date" sidebar, so all the
content could be accessed chronologically.
I understand their comments, though I'm not convinced that a
single, monolithic page of information is the right direction either.
<p>
One feature of blogger that proved helpful was its "Pages", static
html files or links that appear (by default) at the top of the blog.
These links are always present no matter where in the blog you navigate
to. We added the syllabus and course calendar here. And here's where
things start to get interesting.
<p>
All four instructors need the ability to collaboratively edit and then publish
public-facing material with students (like the syllabus), as edit private course material,
like grades. Dropbox makes a natural choice for private material,
whereas a collaborative source repository like GitHub makes a natural
choice for public material. Personally, I'm more familiar/comfortable
with <a href="http://code.google.com/hosting/">Google Code</a>, but I
don't think the choice of services here is material. [Sidenote:
If you have a pro GitHub account, the choice seems easy: use a private GitHub repo.]
<p>
Using a public git repo gave us a few things that I really
liked. First off, our class schedule was a simple HTML table in the
repo that we pointed at with a static link via the Pages widget, above.
When I needed to update an assignment due date, I did a quick
edit-and-commit on my laptop, pushed it to the archive, and magically appeared
on the blog.
If several instructors are simultaneously working on course material,
this essentially provides an audit trail of who did what. This is kindergarden-level
stuff for software developmers. But these are tools that teachers could benefit from
and aren't very familiar with. For example, getting
colleagues to set up git represents a non-trivial challenge. Nonetheless,
there are other benefits of learning git if you're a scientist (that I won't address
here, though see <a href="http://helmingstay.blogspot.com/2013/08/everyday-revision-control.html">here</a> for some thoughts).
<p>
In our class, I also used R for data visualizations.
I let the awesome <a
href="http://yihui.name/knitr/">knitr</a> package build the
appropriate .png files from my R code (along with pdfs for download,
if needed). Again, it sounds simple, but adding the generated figures
to the git archive allowed me to quickly link to them in blog posts,
and then update them later if needed.
<p>
I would have preferred having the blog posts themselves under revision
control (only "Pages" can point to an
external source html without javascript, which I didn't have time for).
But, for the simplicity of setup and the low
(e.g. free) cost of use, I didn't find posts to be much of an issue.
Blogger allows composition in pure html, without all the *junk*, which helps.
But having to re-upload and link a figure every time I find an error?
Definitely a pain.
<p>
For this class, we also set up an <a
href="http://www.e4ward.com/">e4ward.com</a> address that pointed to my
email box. This allowed me to publicly post the class address without
fear of being spammed forevermore, and allowed me to easily identify
all class email. Having a single instructor responsible for email
correspondence worked well enough. As an early assignment, students
were asked to send a question to the class email address. This
is a nice way to get to know folks, and incentivize them to go to
"digital office hours", e.g. ask good questions over email.
We did have some issues with <a
href="http://www.e4ward.com/">e4ward.com</a> towards the end of the class.
This *sucks* - students panic when emails get lost, and it's hard to
sort out where the problem is. Honestly, I don't know the answer here.
<p>
As a sidenote, I pushed the use of Wikipedia (WP) heavily in this
class, and referred to it often myself. This is not possible in all
fields, but WP articles in many of the hard sciences are now both
technically detailed and accessible. Probability and statistics
articles are some of the best examples here, since the "introductory
concept" articles are used by a large number of individuals/fields,
and aren't the subjects of debate (compare this with the WP pages of,
for example, Tibet or Obama). I also discovered WP's "Outlines"
pages for the first time. If you ever use statistics, I strongly
recommend spending some time with WP's <a
href="http://en.wikipedia.org/wiki/Outline_of_statistics">Outline of
Statistics</a>. It's epic.
<p>
One final tool we used quite a bit was <a
href="http://www.limesurvey.org/en/">LimeSurvey</a>. As I was
planning the class, I went looking for an inexpensive yet flexible
survey platform. I was generally disappointed by the offerings. My requirements included
raw data export and more than 10 participants; these features tend
*not* to be free, and survey tools can get pricey. Enter
LimeSurvey (henceforth LS). It's open source, well-documented,
versatile, and simple to use. I was reluctant to invest *too* much
effort in tools I wasn't sure we'd use, but I got LS running in less
than 2 hours. To be fair, our lab has an on-campus, always-on server,
and apache2 is already installed and configured. This would have been
an annoying out-of-pocket expense had I needed to rent a VPS, though
you can now get a lot of VPS for $5/mo. [sidenote: getting our campus
network to play nice with custom DNS was a whole other issue...]
<p>
LS allowed me to easily construct data entry surveys, allowing each
student to enter, for example, their 25 flips of 3 coins, or their
sequence of wins and losses playing the Monty Hall problem. Students
can quickly enter this data from their browser at their convenience.
At its best, visualizations of the class data can give students a sense
of ownership and purpose of in-class exercises.
LS also allowed us to conduct initial background surveys, as well as
anonymous, customized class "satisfaction" surveys mid-course to find
out what was working and what wasn't. We ran into a few administrative
issues with LS, but it's an overall powerful and stable data
collection platform. Students seemed happy with it, and it
provided us with valuable feedback.
<p>
What would I do differently? I failed to set up an email list early
on. It would have been trivial using, e.g. <a
href="https://groups.google.com/">Google Groups</a>. At the time, it
seemed redundant, both for the students and us. Didn't they already
have the blog? In retrospect, it would have proved useful at several
points to communicate "New Blog Post" or "Due Date Changed, check
calendar". I've learned this semester to expect that spoken instructions
are not necessarily heeded, and that receiving
administrative details in writing from multiple sources (blog,
calendar, mailing list) is a Good Idea ™.
<p>
Along this vein, a minor needed improvement is breaking making a
separate table for assignments. I originally combined
assignments with the course calendar in the interest of simplicity, giving
students all the relevant information in one place. In the end, it
was just confusing.
<p>
Overall, the class went very well. We have received positive feedback
from the students, and we now have a detailed digital record of the
course. PDFs of their final project posters are now in the
archive, <a href="https://probability-for-scientists.googlecode.com/git/posters/">where they will live in perpetuity</a>. Personally, I'm not
ready to teach a <a
href="http://en.wikipedia.org/wiki/Massive_open_online_course">MOOC</a>
yet, but I'm sold on digital tools as useful supplements to in-class
material. They allowed me to spend less time doing more. These tools helped
the class "organize itself", and reduced communication overhead between instructors.
To older or less technically-inclined teachers, some of this
might seem difficult or confusing. On one hand, you really only need *one* instructor
to coordinate the tech - it's not very hard for instructors to use the tools once they're set up.
On the other hand, some of the above
pieces are very easy to implement, and support from a department
administrator or technically-inclined teaching assistant (TA) might be
available for more challenging pieces. An installation of LS shared
across a department, for example, should be trivial for any department
that has its own linux server (e.g. math, physics, chem).
<p>
In my mind, a key goal here is that technology not get in the way.
Many of the commercial "online learning products" that I've seen
adopted by universities make simple task complicated, while lacking the flexibility.
In trying to do everything
for everyone, they often end up doing nothing for anyone (or take a
high level of skill or experience to use effectively). I far
preferred using several discrete tools, each of which does a single
job well (class email, blog to communicate, git repo to hold files).
<p>
Are there any interesting tools worth trying out next time? I wonder
how a class twitter-hashtag would work...xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com7The University of New Mexico, Albuquerque, NM 87131, USA35.0843441 -106.6197553000000135.058359100000004 -106.66009580000001 35.1103291 -106.57941480000001tag:blogger.com,1999:blog-4695461192314112143.post-37192562019721310412013-11-26T19:50:00.001-07:002013-12-10T19:08:36.878-07:00The Art of the AlbumWhen I'm working, I'll easily listen to 6 or more hours of music a day. Ten years ago, I listened to a *lot* of public radio, which broadened my ear a lot and introducing me to genres ranging from classical jazz to Native American and <a href="http://kanwstore.myshopify.com/">New Mexican music</a>. Internet radio gave me more choice of stations (I'm now a happy KEXP micro-donor). Finally, there was Rdio (or Spotify, take your pick). Residents of the U.S. got to hear about the wonders of these all-you-can-eat music-as-a-utility services years before they became available here. Complex licensing agreements had to be signed with our musical overlords. But the future has finally arrived, and for $10/month I now have an internet full music (including offline access on phone, $5/mo for just wired computer). <br />
<br />
The magic of a high-quality, easily searched and streamed music archive has transformed the way I listen to music. When I hear a song I like, it now takes me less than 30 seconds to find the album and begin playing on my office computer or phone.<br />
There are a few drawbacks - not every album or artist is available (Joanna Newsom is a particularly galling example), and occasionally I find myself without a reliable cellphone or WiFi signal. But these are minor issues. Overall, the ultimate convenience, the *modern-ness* of it all still blows my mind. To me, this is better than a flying car (of course, I don't even own a normal car). And this convenience has, in the last few years, rekindled my love of the art of the album.<br />
<br />
It seems to me that independent music in general has benefited from digital distribution by allowing artists to more easily break from the more conventional constraints of genre. I see a lot of experimentation here, running all the way up to the Dirty Projectors' avant-garde classical composition style. Growing up in the 90's, I enjoyed Pearl Jam and Nirvana well enough, but much of the "alternative" music that I listened to at the time sounded (and still sounds) rather similar to my ears, e.g. Grunge. The ones that sounded different really stood out, and I still cherish them for it (I'm looking at you, Pixies). Maybe I'm biased now by access to more music and better DJs, but I find the modern American music scene incredibly vibrant and diverse. Every month, I can look forward to new releases from favorite artists, as well as finding something or someone new to open my eyes and make my day.<br />
<br />
What follows is an unordered list of albums that I've recently developed a strong relationship with. These albums cover a wide range of the acoustic/electronic spectrum. I enjoy repetitive, energetic music when I'm working or juggling or cleaning; I love the emotion and classic song-writing of "folk" and "country"; and I love the driving anthems of modern indie. Consequently, I like to think there's "something for everyone" here. And each of these is an *album*, a free-standing work of art worthy of repeated enjoyment in its uncensored, unedited entirety.<br />
<br />
<h1>Obvious:</h1><h3>Macklemore & Ryan Lewis - The Heist (2012)</h3>I'd like to find more music like this: the crossroads between pop and hip hop, independent music that gets radio play, catchy but meaningful. I can think of half a dozen song lyric lines that make great life slogans. Tis is a great album to blast in the car on a warm spring day.<br />
<br />
<h3>alt-J - An Awesome Wave (2012)</h3>I think of alt-J as the Neutral Milk Hotel of this decade: where the hell did they come from? It's such a beautiful, subtle album that came out of nowhere and bears repetition very well. I beg the gods for more in the future.<br />
<br />
<h3>Daft Punk - Random Access Memories (2012)</h3>I never really got into Daft Punk before this album, and I didn't even like it that much the first few plays. The songs on this album tend towards longish, some of them are slow, and I found myself getting bored. Then I began getting lines stuck in my head, and began dipping back in. In the end, I find this an immensely satisfying sort of pop-EDM-concept-album: a soothing mix of repetitive riffs that aren't too fast or insistent with a backdrop of pop anthem melodies. It strikes me as easy-listening Moby? This album is a little slow form me to "sit down" and listen to, but I find it excellent clean-the-house/driving music.<br />
<br />
<h3>Phosphorescent - Muchacho (2013)</h3>Rainy day + hot coffee. Sunset and a beer. Just got dumped, fired, graduated, engaged? This is such an extraordinarily luscious, eloquent album. It makes me remember that I have emotions. Lots of them.<br />
<br />
<h3>Santigold - Master of My Make-Believe (2012)</h3>I always perk up when I hear Santigold singles on the radio, but I was slow to listen to the albums. I like her self-title 2008 album, but it never really got under my skin. The second or third listen of Master, though, and I wanted to know more about this artist. After digging around a bit, I feel like I have a better idea of where she's coming from, and where she's going. The comparison to M.I.A. is inevitable, while the album art for Master suggests something more like Outkast. Master has tons of energy and is packed with pop-friendly riffs. But it's complex, and strikes me as walking the "don't define me" tightrope (or slackline, if you will; you can push *back* on a slackline). I enjoy that it doesn't settle down into a niche and stay there.<br />
<br />
<h1>Less Obvious:</h1><h3>Shovels and Rope - O' Be Joyful (2012)</h3>In my mental map of Americana, I file this near Wilco and Drive By Truckers. Sometimes slow and sweet, sometimes fast and rambunctious, but always melodic, this album is full of luscious 2-part harmonies with a low-fi, intimate feel. I'm always sad when it ends; I always want more.<br />
<br />
<h3>First Aid Kit - The Lion's Roar (2012)<br />
</h3>Can I call this indie-Americana? Less of the overt Southern influence of Shovels and Rope, but still full of tight vocal harmonies of country/folk. Apparently they're sisters, and apparently they're young, but this album has a big sound, full of driving melancholy. Playing two or three of their albums back-to-back is particularly satisfying. They seem to be growing as they go, and I'm excited to hear what comes next.<br />
<br />
<h3>John Grant - Queen of Denmark (2010)<br />
</h3>A very good anthem album. I don't often listen or pay attention to lyrics, but Grant has a John Prine-ish storytelling quality, a dark sense of humor and playful irony. Musically, it's tends towards simple, with a fast, light quality that reminds me of Paul Simon's Graceland. Thematically, though, it's a dark album. A far-off hint of redemption shines at the end of the tunnel, but just barely. Whistling in the dark.<br />
<br />
<br />
<h3>Sharon Van Etten - Tramp (2013)<br />
</h3>A powerful voice, and a powerful song-writer. This album is mature and intimate, and Van Etten's voice is strong and clear. Tight harmonies and vocal stylings that are luxurious without being excessive. The utterly enrapturing quality of controlled liquid of her voice reminds me a little of the Cowboy Junkies' Margo Timmins, with a bit of Joni Mitchell. In short, she's good.<br />
<br />
<h3>Matthew Dear - Beams (2012)<br />
</h3>My first introduction to Matthew Dear, this album is driving. Repetitive, almost grinding, the samples remind me of smoothed-out, slowed down industrial, or gears-and-grease voodoo. It reminds me of being in the belly of a very large machine. The tone palette is less pure than, say, Daft Punk, with lots of glitches and grinding noises. It's also harmonic, full of discordant melodies. And I *love* it. There are songs that I would love to hear on a dark dance floor in a small, crowded night-club. It's sexy as hell, with a floating touch of loss and nostalgia.<br />
<br />
<h3>Jagwar Ma - Howlin (2013)<br />
</h3>This is a somewhat confusing album. A mix of upbeat chorus-driven pop tunes and beat-and-sample driven pop-EDM, I find it a little schizophrenic at times. In the space of two songs, it goes from an drivingly upbeat guitar-and-vocals sound akin to Django Django's recent album, to something more akin to Caribou's hypnotic samples, with little in the way of transition. The situation reminds me a little of Hot Chip's recent album In Our Heads (which I still find deeply confusing). But Howlin is infectious throughout, with several singles that belong in the "party mix". <br />
<br />
<br />
<h3>Dirty Projectors - Swing Lo Megellan (2012)<br />
</h3>Beautiful melodies with a glorious sheen of tightly-controlled noise and discord, this approaches classical composition in broad-scale interest and ability to scare off pedestrians at a first listen. There's just enough rhythm and melody, though, to reel a music-lover in until one gains some familiarity with the subtleties. Then the album really starts to open up. To my ear, it's the opposite of a show-stopping dramatic pop album. It's playful and light, and strange, and curious and coy, going from simple to huge and back. It's complex and, sometimes subtly, very satisfying. This is real sit-down-and-listen music, kind of like going to see the symphony.<br />
<br />
<h3>Junip - Self Title (2013)<br />
</h3>It's not unlikely that you've already heard "Your Life Your Call". I'm sure it's in some movie or another, or will be soon. I get shivers every time I hear this song - like the soundtracks of the Breakfast Club and Trainspotting had a mutant child. Jose Gonzales has a number of solo albums (I'm quite fond of his 2005 Veneer - see below), though I never made the connection with Junip myself. His voice is as clear and emotional as ever, but the sound is bigger and more nuanced, a wonderful blend of semi-acoustic and smooth electronic sounds. This is an emotional album - not any *particular* emotion, but all of them, simultaneously, and a lot. Much like Muchacho, listening to it makes me feel decidedly and acutely human.<br />
<br />
<br />
<h3>Yppah - Eighty One (2012)</h3>Driving indie dream-pop, Yppah's sound is reminiscent of Heliosequence with drum machines. Something to get the shoe-gazers moving around!<br />
<br />
<h1>Less new<br />
</h1>but recently discovered or especially noteworthy, albums follow.<br />
I'm ready to wrap this post up, so these get just a brief mention, but they're all worth a good listen.<br />
<br />
<h3>Caribou - Swim (2010)</h3>Smooth, fast, steady electronic. A masterful album.<br />
<br />
<h3>Jose Gonzales - Veneer (2005)</h3>Contains a cover of The Knife's song "Heartbeats" that I adore. Close and intimate and lush.<br />
<br />
<h3>Crystal Castles - Self Title (2008)</h3>One of my current favorite albums. I think of it as glitch-rock. It's more syth-y than punk, but has a lot of similar aesthetic sensibilities: loud, abrasive, driving, and inspiring. I particularly like to cue up all 3 Crystal Castles albums and listen to them in a go. Loudly.<br />
<br />
<h3>Gold Panda - Lucky Shiner (2010)</h3>Very smooth, incredibly-produced electronic music. Deeply satisfying, good work music. <br />
<br />
<h3>Franz Ferdinand - Self Title (2004)</h3>Anthemic indie-pop. I'm familiar with most of these songs, and was amazed that that they all came from a single album. Buddy Holly meets Lou Reed?<br />
<br />
<h3>Jolie Holland - Escondida (2004)</h3>Lead singer of The Be Good Tanyas, Jolie Holland's solo album is intimate and enwrapping. <br />
<br />
<h3>Juana Molina - Tres Cosas (2004)</h3>Quiet and playful yet insistent percussion is the constant backdrop against which Molina's voice plays. And is it ever playful. Her unassuming Spanish is hypnotizing. She has a new release out that I haven't digested yet, but here's another case where I happily queue up 2 or 3 albums in a row and let them blend effortlessly from one to the next.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-79490359294203227482013-08-23T19:37:00.000-07:002013-08-23T20:10:20.287-07:00Faster, better, cheaper: what is the true value of a computer? One thing I've had a lot of time to think about over the last 15 years is what exactly does faster & more powerful mean? After a decade of clockspeed wars, we've moved on to more cores, more RAM, longer battery life, less weight, backlit keyboards, etc. A new computer still costs about the same as it did 15 years ago... is it any better?<br />
<br />
The more time I spend with the machines, the more I think about usability. A machine is only as powerful as the tasks it can accomplish. I have a $200 netbook (w/linux, of course) that excels at being light. It works great for travel but not netflix. I could not write a paper on it, but I can check email, upload photos/files, etc.<br />
<br />
In our computer cluster here, we upgrade as we hit limits. Ran out of RAM? Buy more... it's cheap (for a desktop, at least). Chronic, awful wrist pain? Get an ergonomic keyboard. I find a 2-screen setup a very cost-effective productivity boost, whereas the idea of paying 20% more for a 10% bump in speed strikes me as silly.<br />
<br />
The whole state of affairs reminds me a little of mean and variance. We always hear about the expected values of things, the mean, but rarely is variance ever reported, and variance is often the most important part. Things like weather, lifespan, salary, time-to-completion? The variance may be more important than the mean. Speed/power tells something about the "maximum potential" of a machine, but not how much use one might get out of it.<br />
<br />
I see two different directions --<br />
<br />
First, unexpected developments. Examples include multiple cores and SSD. No one really expected these to change the aesthetics of computers. Nonetheless, the reduced latency of SSD is pleasantly surprising, as is the increased responsiveness under load of a multi-core machine. I don't hit enter and wait. My machine flows more at the speed of thought than it used to, even if I have a browser with 50+ tabs open, music playing, etc.<br />
<br />
Second, human interface. My new phone is just a little too tall, which makes it just a little more difficult to use, since I can't reach across the screen with one hand. Again, I never expected this to matter, but it's a human-machine interface question rather than a pure machine capabilities question. Backlit keyboards, ergonomic keyboards, featherweight laptops, and long battery lives are all about the human-machine interface. Which is more important than speed, since the human is the whole point....<br />
<br />
The aesthetics of interface is how Apple came to rule the world. Their hardware is beautiful, intuitively responsive to human touch. Hell, even their stores are clean & informative and intuitive and full of fun toys. Personally, I can't stand their walled-garden approach to hardware, and I have the time and energy to coax my machines into greatness through a commandline interface (which remains one of the most powerful human-machine interface ever developed), but I *like* Apple hardware. They've dragged the PC industry kicking and screaming into the 21st century of "Humans matter more than machines".<br />
<br />
Which they rather presciently highlighted in their 1984 Superbowl commercial:<br />
<a href="http://www.youtube.com/watch?v=2zfqw8nhUwA">http://www.youtube.com/watch?v=2zfqw8nhUwA</a><br />
<br />
Finally, a mind-blowing historical view on the subject, including a 1983 Bill Gates pimping Apple hardware, and Steve Jobs describing how machines should help people rather than the other way around:<br />
<a href="http://www.youtube.com/watch?v=xItV5U-V2W4">http://www.youtube.com/watch?v=xItV5U-V2W4</a>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-68668546925920604122013-08-23T19:12:00.000-07:002013-08-23T19:12:02.047-07:00Everyday revision control This post has been a long time coming. Over the past year or so, I've gradually become familiar, even comfortable, with <a href="http://en.wikipedia.org/wiki/Git_%28software%29">git</a>. I've mainly used it for my own work, rather than as a collaborative tool. Most of the folks that I work with don't need to share code on a day-to-day basis, and there's a learning curve that few of my current colleagues seem interested in climbing at this point. This hasn't stopped me from *talking* to my colleagues about git as an important tool in reproducible research (henceforth referred to as RR). <br />
<br />
I find the process of committing files and writing commit messages at the end of the day forces me to tidy up. It also allows me to more easily put a project on hold for weeks or months and to then return to it with a clear understanding of what I'd been working on, and what work remained. In short, I use my git commit messages very much like a lab notebook (a countervailing view on git and reproducible research is <a href="http://kbroman.wordpress.com/tag/reproducible-research/">here</a>, an interesting discussion of GNU make and RR <a href="http://kbroman.github.io/minimal_make/">here</a>, and a nice post on RR from knitr author Yihui Xie <a href="http://www.r-bloggers.com/making-reproducible-research-enjoyable/">here</a> ).<br />
<br />
<br />
Sidenote: I've hosted several projects at <a href="https://code.google.com/">https://code.google.com/</a>, and used their git archives, particularly for classes (I prefer the interface to github, though the two platforms are similar). I've also increasingly used Dropbox for collaborations, and I've struggled to integrate Dropbox and Git. Placing the same file under the control of two very different synchronization tools strikes me as a Bad Idea (TM), and Dropbox's handling of symlinks isn't very sophisticated. On the other hand, maintaining 2 different file-trees for one project is frustrating. I haven't found a good solution for this yet...<br />
<br />
As far as tools go, most of the time I simply edit files with vim and commit from the commandline. In this sense, git has barely changed my work flow, other than demanding a bit of much-needed attention to organization on my part. Lately, I've started using GUI tools to easily visualize repositories, e.g. to simultaneously see a commit's date, message, files, and diff. Both gitk and giggle have similar functionality -- giggle is prettier, gitk is cleaner. Another interesting development is that Rstudio now includes git tools (as well as native latex and knitr support in the native Rstudio editor). This means that a default Rstudio install has all the tools necessary for a collaborator to quickly and easily check out an repository and start playing with it.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-23194968509014775872013-08-09T04:20:00.001-07:002013-08-09T04:25:10.902-07:00Adventures with AndroidAfter months of dealing with an increasingly sluggish and downright buggy Verizon HTC Rhyme, I finally took the leap and got a used Galaxy Nexus. First off, I think it's beautiful. The rhyme isn't exactly a high-end phone, so the small, unexciting screen isn't particularly surprising. By comparison, the circa 2011 Nexus is a work of art. My first impression of the AMOLED screen is great. Dark blacks and luscious color saturation (though I have found it to be annoyingly shiny -- screens shouldn't be mirrors!). It even has a <a href="http://gizmodo.com/5851288/why-the-barometer-is-androids-new-trump-card">barometer</a>!<br />
<br />
One big motivation for a phone upgrade (asides from the cracked screen and aforementioned lag) was being stuck at Android 2.3.4 (Gingerbread). As a technologist, I don't consider myself an early-adopted. I prefer to let others sort out the confusion of initial releases, and pick out the gems that emerge. But Gingerbread is well over 2 years old, and a lot has happened since then. The Galaxy Nexus (I have the Verizon CDMA model, codename <a href="http://wiki.cyanogenmod.org/w/Toro_Info">Toro</a>) is a skin-free, pure Android device. Which means that I am now in control of my phone's Android destiny!<br />
<br />
How to go about this? I hit the web and cobbled together a cursory understanding of the Android/google-phone developer ecosystem as it currently stands. First off, there's <a href="http://forum.xda-developers.com/">xda-developers</a>, a very active community of devs and users. There's an organizational page for information on the Galaxy Nexus <a href="http://forum.xda-developers.com/showpost.php?p=32317030&postcount=2">here</a> that helped me get oriented. <a href="http://www.webupd8.org/2012/08/install-adb-and-fastboot-android-tools.html">This post</a> made installing adb and fastboot a snap on ubuntu 12.04 (precise). There's also some udev magic from google <a href="http://source.android.com/source/initializing.html"> under the "Configuring USB Access" section here</a> that I followed, perhaps blindly (though http://source.android.com is a good primary reference...).<br />
<br />
Next, I <a href="http://clockworkmod.com/rommanager">downloaded ClockworkMod</a> for my device, rebooted my phone into the bootloader, and installed and booted into ClockworkMod:<br />
<br />
<pre>## these commands are from computer connected to phone (via usb cable)
## check that phone is connected. this should return "device"
adb get-state
## reboot the phone into the bootloader
adb reboot-bootloader
## in recovery, the phone should have funny pictures of the adroid bot opened up... reminds me of Bender.
## the actual file (the last argument) will vary by device
fastboot flash recovery recovery-clockwork-touch-6.0.3.5-toro.img
## boot into ClockworkMod
fastboot boot recovery-clockwork-touch-6.0.3.5-toro.img
</pre><br />
This brought me to the touchscreen interface of ClockworkMod. First, I did a factory reset/clear cache as per others' instructions. Then I flashed the files listed <a href="http://forum.xda-developers.com/showthread.php?t=2395682">here</a> (with the exception of root) via sideloading. There's an option in ClockworkMod that says something like "install from sideload". Selecting this gives instructions -- basically, use adb to send the files, and then ClockworkMod takes care of the rest:<br />
<br />
<pre>## do this on computer after selecting "install from sideload" on phone
## the ROM, "pure_toro-v3-eng-4.3_r1.zip", varies by device
adb sideload pure_toro-v3-eng-4.3_r1.zip
## repeat for all files that need to be flashed
</pre><br />
I rebooted into a shiny new install of Jelly Bean (4.3). It's so much cleaner and more pleasant than my old phone. I was also pleasantly surprised to see that Android Backup service auto-installed all my apps from the Rhyme.<br />
<br />
In the process of researching this, I got a much better idea of what <a href="http://en.wikipedia.org/wiki/CyanogenMod">CyanogenMod</a> does. I'm tempted to try it out now, but I reckon I'll wait for the 4.3 release, whenever that happens.<br />
<br />
<br />
I also found <a href="http://www.htcdev.com/bootloader">http://www.htcdev.com/bootloader</a>, which offers the prospect of unlocking and upgrading the HTC Rhyme, though I haven't found any ROMs that work for the CDMA version... xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-61587104934733115792013-06-13T04:18:00.001-07:002013-06-13T04:18:56.270-07:00Secure webserver on the cheap: free SSL certificatesSetting up an honest, fully-certified secure web server (e.g. https) on the cheap can be tricky, mainly due to certificates. Certificates are only issued to folks who can prove they are who they say they are. This verification generally takes time and energy, and therefore money. But the great folks at <a href="https://www.startssl.com/">https://www.startssl.com/</a> have an automated system that verifies identity and auto-renders associated SSL certificates for free.<br />
<br />
Validating an email is easy enough, but validating a domain is trickier -- it requires a receiving mailserver that startssl can mail a verification code to. Inbound port 25 (mail server) is blocked by <a href="http://mailreg.unm.edu/">my ISP, the University of New Mexico</a> (and honestly, I'd rather not run an inbound mail server).<br />
<br />
I manage my personal domain through <a href="http://freedns.afraid.org/">http://freedns.afraid.org/</a>. They provide full DNS management, as well as some great dynamic DNS tools. They're wonderful. But they don't provide any fine-grained email management, just MX records and the like.<br />
<br />
The perfect companion of afraid.org is <a href="https://www.e4ward.com/">https://www.e4ward.com/</a>. They have mail servers that will conditionally accept mail for specific addresses at personal domain, and forward that mail to an email account. This lets me route specific addresses @mydomain.com, things like postmaster@mydomain.com, to my personal gmail account. E4ward is <a href="http://newhelp.e4ward.com/">a real class-act</a>. They manually moderate/approve new accounts, so there's a bit of time lag. To add a domain, they also require proof of control via a TXT record (done through afraid.org).<br />
<br />
This whole setup allowed me to prove that I owned my domain to startssl.com without running a mail server or paying for anything other than the domain. The result is my own SSL certificates. I'm running a pylons webapp with apache2 and mod_wsgi. In combination with python's repoze.what, I get secure user authentication over https without any <a href="http://askubuntu.com/questions/49196/how-do-i-create-a-self-signed-ssl-certificate">snakeoil</a>.<br />
<br />
Hat-tip to <a href="http://forum.asb-sakray.net/index.php?showtopic=32242">this writeup</a>, which introduced me to e4ward.com and their mail servers.<br />
<br />
Finally, there are a number of online tools to query domains. <a href="http://www.dnsstuff.com/tools#dnsReport">dnsstuff.com</a> was one of the better ones I found. It takes longer to load, but gives a detailed report of domain configuration, along with suggestions. A nice tool to verify that everything is working as expected. <br />
<br />
<br />
xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com5tag:blogger.com,1999:blog-4695461192314112143.post-21559408497976437312013-06-11T04:27:00.000-07:002013-06-11T04:37:23.609-07:00Learning new fileserver tricks: RAID + LVMI've finally <a href="http://www.ducea.com/2009/03/08/mdadm-cheat-sheet/">gotten comfortable with linux's software raid</a>, aka mdadm. I've been hearing about LVM, and I finally took the plunge and figured out how to get the two to play together. Of course, a benefit of RAID is data security. The big benefit I see from LVM is getting to add/remove disk space without repartitioning. Once RAID is working, stacking LVM on top was easy enough, especially for my use case of a single-big-filesystem. I was able to move all my data onto one RAID array, built a new filesystem on top of a logical volume, move data to the new filesystem, and then add the final RAID array to the logical volume and resize the filesystem. Thus, I end up with 3 separate RAID arrays glommed together into a single, large filesystem. <br />
<br />
<pre>## Tell LVM about RAID arrays
sudo pvcreate /dev/md2
sudo pvcreate /dev/md3
## Create a volume group from empty RAID arrays
sudo vgcreate VolGroupArray /dev/md2 /dev/md3
## Create a logical volume named "archive", using all available space
sudo lvcreate -l +100%FREE VolGroupArray -n archive
sudo lvdisplay
## and create a filesystem on the new logical volume
sudo mkfs.ext4 /dev/VolGroupArray/archive
## mount the new filesystem
## and move files from the mount-point of /dev/md1 to /dev/VolGroupArray/archive
## then unmount /dev/md1
## Add the last RAID array to the volume group
sudo pvcreate /dev/md1
sudo vgextend VolGroupArray /dev/md1
## Update the logical volume to use all available space
sudo lvresize -l +100%FREE /dev/VolGroupArray/archive
## And resize the filesystem -- rather slow, maybe faster to unmount it first...
sudo resize2fs /dev/VolGroupArray/archive
## Finally, get blkid and update /etc/fstab with UUID and mount options (here, just noatime)
sudo blkid
</pre><br />
I probably should have made backups before I did this, but everything went smoothly...<br />
Also, I discovered <a href="https://github.com/g2p/blocks">this python tool</a> to do conversions in-place. Again, this appears non-destructive, but back-ups never hurt. Also of interest for a file server is Smartmontools to monitor for hardware/disk failures: a nice review <a href="https://help.ubuntu.com/community/Smartmontools"> is here</a>.<br />
<br />
[REFS]<br />
* <a href="http://home.gagme.com/greg/linux/raid-lvm.php">http://home.gagme.com/greg/linux/raid-lvm.php</a><br />
* <a href="https://wiki.archlinux.org/index.php/Software_RAID_and_LVM">https://wiki.archlinux.org/index.php/Software_RAID_and_LVM</a><br />
* <a href="http://webworxshop.com/2009/10/10/online-filesystem-resizing-with-lvm">http://webworxshop.com/2009/10/10/online-filesystem-resizing-with-lvm</a>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-43620929023508890862013-06-07T01:28:00.001-07:002013-06-07T01:37:40.267-07:00Symmetric set differences in RMy .Rprofile contains a collection of convenience functions and function abbreviations. These are either functions I use dozens of times a day and prefer not to type in full:<br />
<pre>## my abbreviation of head()
h <- function(x, n=10) head(x, n)
## and summary()
ss <- summary
</pre>
Or problems that I'd rather figure out once, and only once:
<pre>## example:
## between( 1:10, 5.5, 6.5 )
between <- function(x, low, high, ineq=F) {
## like SQL between, return logical index
if (ineq) {
x >= low & x <= high
} else {
x > low & x < high
}
}
</pre>
One of these "problems" that's been rattling around in my head is the fact that setdiff(x, y) is asymmetric, and has no options to modify this.
With some regularity, I want to know if two sets are equal, and if not, what are the differing elements. setequal(x, y) gives me a boolean answer to the first question. It would *seem* that setdiff(x, y) would identify those elements. However, I find the following result rather counter-intuitive:
<pre>> setdiff(1:5, 1:6)
integer(0)
</pre>I personally dislike having to type both setdiff(x,y) and setdiff(y,x) to identify the differing elements, as well as remember which is the reference set (here, the second argument, which I find personally counterintuitive). With this in mind, here's a snappy little function that returns the symmetric set difference:
<pre>symdiff <- function( x, y) { setdiff( union(x, y), intersect(x, y))}
> symdiff(1:5, 1:6) == symdiff(1:6, 1:5)
[1] TRUE
</pre><br />
Tada! A new function for my .Rprofile!xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com1tag:blogger.com,1999:blog-4695461192314112143.post-92196005097788352162012-11-25T00:44:00.000-07:002012-11-25T01:16:18.805-07:00Successive Differences of a Randomly-Generated TimeseriesI was wondering about the <i>null distribution</i> of successive differences of random sequences, and decided to do some numerical experiments. I quickly realized that successive differences equates to taking successively higher-order numerical derivatives, which functions as a high-pass filter. So, the <i>null distribution</i> really depends on the spectrum of the original timeseries.
<br>
<br>
Here, I've only played with random sequences, which are equivalent to white noise. Using the wonderful animation package, I created a movie that shows the timeseries resulting from differencing, along with their associated power spectra. You can see that, by the end, almost all of the power is concentrated in the highest frequencies. The code required to reproduce this video is shown below.
<br><br>
Note -- For optimum viewing, switch the player quality to high.
<iframe width="800" height="600" src="http://www.youtube.com/embed/y9jKJnF0IJI?hd=1" frameborder="0" allowfullscreen></iframe>
<pre>
require(animation)
## large canvas, write output to this directory, 1 second per frame
ani.options(ani.width=1280, ani.height=720, loop=F, title='Successive differences of white noise', outdir='.', interval=1)
require(plyr)
require(xts)
## How many realizations to plot?
N=5
## random numbers
aa = sapply(1:26, function(x) rnorm(1e2));
colnames(aa) = LETTERS;
saveVideo( {
## for successive differences, do...
for (ii in 1:50) {
## first make the differences and normalize
aa1 = apply(aa, 2, function(x) {
ret=diff(x, differences=ii);ret=ret/max(ret)
});
## Turn into timeseries object for easy plotting
aa1 = xts(aa1, as.Date(1:nrow(aa1)));
## next, compute spectrum
aa2 = alply(aa1, 2, function(x) {
## of each column, don't modify original data
ret=spectrum(x, taper=0, fast=F, detrend=F, plot=F);
## turn into timeseries object
ret= zoo(ret$spec, ret$freq)});
## combine into timeseries matrix
aa2 = do.call(cbind, aa2 );
colnames(aa2) = LETTERS;
## plot of timeseries differences
## manually set limits so plot area is exactly the same between successive figures
myplot = xyplot(aa1[,1:N], layout=c(N,1),
xlab=sprintf('Difference order = %d', ii),
ylab='Normalized Difference',
ylim=c(-1.5, 1.5),
scales=list(alternating=F, x=list(rot=90), y=list(relation='same')));
## plot of spectrum
myplot1 = xyplot(aa2[,1:N], layout=c(N,1),
ylim=c(-0.01, 5), xlim=c(0.1, 0.51),
xlab='Frequency', ylab='Spectral Density',
type=c('h','l'),
scales=list(y=list(relation='same'), alternating=F));
## write them to canvas
plot(myplot, split=c(1,1,1,2), more=T);
plot(myplot1, split=c(1,2,1,2), more=F);
## provide feedback of process
print(ii)}
## controls for ffmpeg
}, other.opts = "-b 5000k -bt 5000k")
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-27720233664554958962012-11-08T23:36:00.000-07:002012-11-08T23:40:03.621-07:00Consumerist DilemasConsumer society gives us lots of choices, but sometimes provides little opportunity for post-choice customization. We buy something, and we can either return it or keep it, but it is what it is. Want something a little shorter, a little tighter, or a little less stiff? Too bad. Unless you're loaded, and then you can hire someone to do a custom job. But, to my knowledge, even the super-rich no longer order custom automobiles.
<br><br>
For a small subset of expensive, long-lived personal items like cars and mattresses, this can make new purchases particularly stressful. This is one plausible explanation for the copious consumer-choice media devoted to, for example, cars, as well as the effort companies put into brand identity. If they sell you something that you like, then there's a reasonable chance you'll get another one, the "safe choice", when the old one breaks.
<br><br>
I've faced this dilemma recently with shoes and glasses. I typically have one or two pairs of each that I use primarily, daily, typically for at least 2 years. I'm very near-sighted, so I have to get the expensive high-index lenses. My big toes spread out, I assume from years of wearing flip-flops and Chacos, so most every shoe feels narrow and pointy. The decision to get a new pair of glasses or shoes fills me with an existential dread -- what if I'm stuck with costly discomfort for the next year or more?
<br><br>
This month I discovered that keyboards engender a similar dread. I'm at a keyboard for at least 20 hours on an average week, and I imagine it gets up to 60 hours on a long week. Not all of that is typing, but I've had a slowly building case of carpal tunnel / repetitive stress injury, particularly in my right hand, from long coding sessions that involve heavy use of shift, up-arrow, and enter keys. After a month of poking around the internet, checking reviews on Amazon and looking at eBay for used items, I felt totally conflicted. I just wanted to "try something on" before I shelled out cash to saddle myself with something that I ended up hating.
<br><br>
But anything would be better than my stock Dell monstrosity.
In this spirit, I dropped by the chaotic indie computer parts/repair store across the street from campus, and found two ergonomic keyboards used and on the cheap:
<br><br> #1 Microsoft Comfort Curve 2000 v1
<br> #2 Adesso EKB-2100W
<br><br>
The perfect keyboard is neither of them. The Adesso has a nice split design, but it's huge. Enormous. Heavy -- it could be an instrument of murder in the game of Clue. #1 looks nice, and has really easy key-travel. Unfortunately, the left Caps-lock, which is mapped to control and is involved in something like ~5-20% of my (non-prose) keystrokes, is sticky and giving my pinky problems already.
<br><br>
As least they're cheap.
<br><br>
A friend pointed out that one key to preventing RSI is to change up the routine. With this in mind, maybe two keyboards isn't such a bad idea. Did I mention they're cheap? At very least, the existential dread has abated. No more keyboard choices to make in the foreseeable future. Now, if I could just find a low-profile pair of sneakers with a large toebox and a tasteful lack of neon and mesh.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com1tag:blogger.com,1999:blog-4695461192314112143.post-86705676370254766402012-10-02T02:56:00.000-07:002012-11-08T23:37:27.421-07:00Modeling PhilosopyI've found it sometimes difficult to explain my modeling philosophy, especially in the face of "why don't you just do X?" or "isn't that just like Y and Z?" responses. Without going into X, Y, or Z in the slightest, I present the following quote from one of my intellectual heros. His short book (less than 100 pages) full of equations that I struggle to make sense of (not unlike going to the opera) is full of little morsels of wisdom. Here's one of my favorites:
<blockquote>It need hardly be repeated that no ~detailed~ statistical agreement with observation is to be expected from models which are based on admitted simplifications. But the situation is quite analogous to others where complex and interacting statistical systems are under consideration. What is to be looked for is a comparable pattern with the main features of the phenomena; when, as we shall see, the model predicts important features which are found to conform to reality, it becomes one worthy of further study and elaboration.</blockquote>
--M.S Bartlett, 1960.
<br>
Stochastic Population Models in Ecology and Epidemiologyxianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-48831072053132543292012-02-29T00:14:00.016-07:002012-03-08T04:17:46.361-07:00Custom Amazon EC2 config for Rstudio<h3>Introduction</h3>This post is a work in progress building on the previous post. It's my attempt to simultaneously learn Amazon's AWS tools and set up R and Rstudio Server on a customized "cloud" instance. I look forward to testing some R jobs that have large memory requirements or are very parallelizable in the future.<br />
<br />
To start, I followed the instructions <a href="http://bioconductor.org/help/bioconductor-cloud-ami/">here</a> to get a vanilla Ubuntu/Bioconductor EC2 image up and running. This was a nice exercise -- props to the bioconductor folks for setting up this image.<br />
<br />
<h4>Edit -- another look </h4>After playing with this setup more and trying to shrink the root partition (as per <br />
<a href="http://www.thatsgeeky.com/2011/11/shrink-ebs-root/">these instructions</a>), I realized there's a tremendous amount of cruft in this AMI. It's over 20GB, there's a weird (big) /downloads/ folder lying aruond, and just about every R package ever in /use/local/lib64/R. I cut it down to ~6gb by removing these two directories. <br />
<br />
<br />
If I were to do this again, I would use the <a href="http://uec-images.ubuntu.com/releases/10.04/release/">official Canonical images</a> (more details <a href="https://help.ubuntu.com/community/EC2StartersGuide">here</a>), and manually install what I need. Honesty, this is more my style -- it means fewer unknown hands have touched my root, and there's a clear audit chain of the image. Aside from installing a few extra packages (<a href="http://rstudio.org/download/server">Rstudio server, for example</a>), the rest of the instructions should be similar.<br />
<br />
<h3>Modifications</h3>I proceeded to lock it down -- add an admin user, prevent root login, etc.<br />
Then I set up apache2 to serve Rstudio server over HTTPS/SSL.<br />
In the AWS console, I edited Security Groups to add a custom Inbound rule for TCP port 443 (https). I then closed off every other port besides 22 (ssh).<br />
<br />
Below is the commandline session that I used to do it, along with annotations and links to some files.<br />
<pre>adduser myusername
adduser myusername sudoers
## add correct key to ~myusername/.ssh/authorized_keys
vi /etc/ssh/sshd_config
## disable root login
/etc/init.d/ssh restart
## now log in as myusername via another terminal to make sure it works, and then log out as root
</pre><br />
<br />
Next, I set up dynamic dns using http://afraid.org (I've previously registered my own domain and added it to their system). I use a script file made specifically to work with AWS -- <a href="http://helmingstay.blogspot.com/p/afraidorg-script-for-amazon-ec2.html">it's very self-explanatory</a>. <br />
<br />
<pre>## change hostname to match afraid.org entry
sudo vi /etc/hostname
sudo /etc/init.d/hostname restart
## Now it's time to make Rstudio server a little more secure
## from http://rstudio.org/docs/server/running_with_proxy
sudo apt-get install apache2 libapache2-mod-proxy-html libxml2-dev
sudo a2enmod proxy
sudo a2enmod proxy_http
## based on instructions from http://beeznest.wordpress.com/2008/04/25/how-to-configure-https-on-apache-2/
openssl req -new -x509 -days 365 -keyout crocus.key -out crocus.crt -nodes -subj \
'/O=UNM Biology Wearing Disease Ecology Lab/OU=Christian Gunning, chief technologist/CN=crocus.x14n.org'
## change permissions, move to default ubuntu locations
sudo chmod 444 crocus.crt
sudo chown root.ssl-cert crocus.crt
sudo mv crocus.crt /usr/share/ca-certificates
sudo chmod 640 crocus.key
sudo chown root.ssl-cert crocus.key
sudo mv crocus.key /etc/ssl/private/
sudo a2enmod rewrite
sudo a2enmod ssl
sudo vi /etc/apache2/sites-enabled/000-default
sudo /etc/init.d/apache2 restart
</pre><br />
You can see my full apache config file <a href="http://helmingstay.blogspot.com/p/apache2-config-for-ssl.html">here</a>:<br />
<br />
<h3>Conclusions</h3>Now I access Rstudio on the EC2 instance with a browser via:<br />
https://myhostname.mydomain.org<br />
<br />
I found that connecting to the Rstudio server web interface gave noticable lag. Most annoyingly, key-presses were missed, meaning that I kept hitting enter on incorrect commands. Connecting to the commandline via SSH worked much better.<br />
<br />
Another annoyance was that Rstudio installs packages into something like ~R/libraries, whereas the commandline R installs them into ~/R/x86_64-pc-linux-gnu-library/2.14. Is this a general feature of Rstudio? It's a little confusing that this isn't standardized.<br />
<br />
Another quirk -- I did all of this on a Spot Price instance. After all of these modifications, I discovered that Spot instances can't be "stopped" (the equivalent of powering down), only terminated (which discards all of the changes). After some looking, I discovered that I could "Create an Image" (EBS AMI) from the running image. This worked well -- I was able to create a new instance from the new AMI that had all of the changes, terminate the original instance, and then stop the new instance.<br />
<br />
All of this sounds awfully complicated. Overall, this is how I've felt about Amazon AWS in general and EC2 in particular for a while. The docs aren't great, the web-based GUI tools are sometimes slow to respond, and the concepts are *new*. But I'm glad I waded in and got my feet wet. I now understand how to power up my customized image on a micro instance for less than $0.10 an hour to configure and re-image it, and how to run that image on an 8 core instance with 50+GB RAM for less than a dollar an hour via Spot Pricing.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-53267364945981792872012-02-28T03:17:00.004-07:002012-02-29T02:53:34.112-07:00Adventures in R Studio Server: Apache2, Https, Security, and Amazon EC2.I just put a fresh install of Ubuntu Server (10.04.4 LTS) on one of our machines. As I was doing some post-install config, I accidentally installed Rstudio Server. And subsequently fell down an exciting little rabbit-hole of server configuration and "ooooh-lala!" playtime.<br />
<br />
A friend sung the wonders of Rstudio Server to me recently, and I filed it under "things to ignore for now". Just another thing to learn, right? Turns out, the Rstudio folks do *great* work and write good docs, so I hardly had to learn anything. I just had to dust off my sysadmin skills and fire up some google.<br />
<br />
I'm a little concerned about running web services on public-facing machines. Even more so, given that R provides fairly low-level access to operating system services. Still, I was impressed to see system user authentication.<br />
<br />
I followed the docs for <a href="http://rstudio.org/docs/server/running_with_proxy">running apache2 as a proxy server</a>, and learned a little about apache in the process. Since I made it this far, I figured I'd run it <a href="http://beeznest.wordpress.com/2008/04/25/how-to-configure-https-on-apache-2/">through https/ssl</a>, add some <a href="http://rstudio.org/docs/server/configuration">memory limitations</a>, etc. I'm still not entirely convinced this is secure -- it seems that running it in a virtual machine or <a href="http://kegel.com/crosstool/current/doc/chroot-login-howto.html">chroot jail</a> would be ideal.<br />
<br />
On the other hand, I ran across <a href="http://toreopsahl.com/2011/10/17/securely-using-r-and-rstudio-on-amazons-ec2/">this post</a> on running Rstudio Server inside Amazon EC2 instances. Nighttime <a href="http://aws.amazon.com/ec2/pricing/">EC2 spot prices </a>on "<a href="http://aws.amazon.com/ec2/instance-types/">Quadruple Extra Large</a>" instances (68.4 GB of memory, <br />
8 virtual cores with 3.25 EC2 Compute Units each) fell below $1 an hour tonight, which is cheap enough to play with for an hour or two -- take it through some paces and see how well it does with a *very* *large* *job* or two. Instances can now be stopped and saved to EBS (elastic block storage), and so only need to be configured once, which really simplifies matters. In fact, I'm wondering if Rstudio (well, R, really) is my "killer app" for EC2. <br />
<br />
Overall, I was really impressed at how fast and easy this was to get up and running. Fun times ahead!xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com1tag:blogger.com,1999:blog-4695461192314112143.post-44049961888709338322011-10-29T00:13:00.002-07:002012-02-29T02:57:01.879-07:00SQL KoanIt's not <b>that</b> profound, but I sure do like it: a simple, elegant example of a self-join that gives a truth table for the <a href="http://wiki.postgresql.org/wiki/Is_distinct_from">NOT DISTINCT FROM</a> operator. <br />
<br />
<pre>WITH x(v) AS (VALUES (1),(2),(NULL))
SELECT l.v, r.v, l.v = r.v AS equality,
l.v IS NOT DISTINCT FROM r.v AS isnotdistinctfrom
FROM x l, x r;
</pre><br />
Thanks to <a href="http://postgresql.1045698.n5.nabble.com/SELF-LEFT-OUTER-JOIN-SELF-JOIN-including-NULL-values-tp2843950p2843982.html">Sam Mason</a> via the PostgreSQL - general mailing list.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-240083545361145112011-10-28T21:42:00.001-07:002012-02-29T02:54:59.597-07:00Why can't ads be more fun?The worst part about advertisements these days, in my not-so-humble opinion, is the degree to which they belittle the user. Well-endowed lady wants to be my Facebook friend? "Click here to confirm request." Riiight. Lonely singles near me? Yeah, sure. <a href="http://www.wired.com/magazine/2011/09/mf_scareware/all/1">OH GOOD LORD, A THREAT HAS BEEN IDENTIFIED!</a> Oh, wait, y'all really <b>are</b> evil! On the infrequent occasions that I do watch TV, the ads are less flagrantly insidious, but they're nonetheless relentlessly patronizing.<br />
<br />
The second-worst part of ads, IMNSHO, is the brute-force repetition. The same inane 30 seconds once? I can do that, but five times in an hour? You've got to be kidding me. Nope, no kidding here. Just 30 seconds of the <b>same</b> inanity, again and again.<br />
<br />
Consequently, I find unique and intelligent ads rather inspiring. I know, it shouldn't be this way. One might thing that unique and inspiring would be more common amongst something so prevalent. Anyway...<br />
<br />
Some of the GEICO ads that have played on <a href="http://hulu.com">Hulu</a> lately have managed to avoid the "worst part". They're short and playful, vaguely reminiscent of <a href="http://en.wikipedia.org/wiki/Adult_swim">[adult swim]</a> commercials. "Is the sword mightier than the pen" is my current favorite; I still sometimes chuckle when I see it. They still fail occasionally in the repetition department. If you're a large international corporation, how hard is it to make more than a handful of moderately-entertaining 30-second spots?<br />
<br />
<br />
Which brings me back to banner ads. I just saw this at the top of <a href="http://arstechnica.com">ars technica</a> today for the first time, and it really caught my eye. I don't read or write <a href="http://en.wikipedia.org/wiki/Javascript">Javascript</a>, but I found myself puzzling through it and then following <a href="http://heinz.cmu.edu/c/heinz-mism-code.html?utm_campaign=mism-june2011&utm_medium=banner&utm_source=wired-ars&utm_content=code">the link</a>. <br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-130qNWKZNiM/TquCFTkJovI/AAAAAAAABC0/UcGtKFRlJSc/s1600/5333546_Ad-2-Banner.jpg" imageanchor="1" style=""><img border="0" height="49" width="400" src="http://4.bp.blogspot.com/-130qNWKZNiM/TquCFTkJovI/AAAAAAAABC0/UcGtKFRlJSc/s400/5333546_Ad-2-Banner.jpg" /></a></div><br />
At the top of the page, I find this banner, which visually reminds me where I came from, that I've come to the right place, and gives me another little puzzle. Neat. Thank you, Heinz College of Carnegie Mellon.<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-3a6i63FsgVQ/TquCL4Rq2xI/AAAAAAAABDA/_7YIUIVuja0/s1600/banner_code.jpg" imageanchor="1" style=""><img border="0" height="59" width="400" src="http://3.bp.blogspot.com/-3a6i63FsgVQ/TquCL4Rq2xI/AAAAAAAABDA/_7YIUIVuja0/s400/banner_code.jpg" /></a></div><br />
None of this is revolutionary, but I <b>do</b> hope it is evolutionary. Ads don't <b>have</b> to be boring, and I'm guessing that interesting and intelligent ads are more effective even as they're less painful and insulting. My guess is that a generational change-over in marketing departments and their managers is underway, and will be slow. Still, I look forwards to a new generation of more modern advertising approaches that gives me something to <b>think</b> while my eyeballs are held hostage.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-77183409073473194972011-10-06T20:28:00.000-07:002011-10-06T20:28:03.207-07:00FFT / Power Spectrum Box-and-Whisker Plot with Gggplot2I have a bunch of time series whose power spectra (FFT via <b>R</b>'s <b>spectrum()</b> function) I've been trying to visualize in an intuitive, aesthetically appealing way. At first, I just used lattice's <b>bwplot</b>, but the spacing of the X-axis here really matters. The spectra's frequencies aren't regularly-spaced categories, which is the default of <b>bwplot</b>. If all the series are of the same length, then the coefficients of the are all estimated at the same frequencies, and one can use <b>plyr</b>'s split-apply-combine magic to compute the distribution of coefficients at each frequency.<br />
<br />
I spent some time trying to faithfully reproduce the plotting layout of the figures produced by, e.g. <b>spectrum(rnorm(1e3))</b> with respect to the axes. This was way more annoying than i expected...<br />
<br />
To include a generic example, I started looking around for easily generated, interesting time series. The <a href="http://en.wikipedia.org/wiki/Logistic_map">logistic map</a> is <a href="http://helmingstay.blogspot.com/2011/10/bat-country.html">so damned easy to compute</a> and generally pretty interesting, capable of producing a range of frequencies. Still, I was rather surprised by the final output. Playing with these parameters produces a range of serious weirdness -- apparent ghosting, interesting asymmetries, etc.. <br />
<br />
Note that I've just used the logistic map for an interesting, self-contained example. You can use the exact same code to generate your own <b>spectrum()</b> box-and-whisker plot by substituting your matrix of time series for <b>mytimeseries</b> below. Do note here that series are by row, which is not exactly standard. For series by column, just change the margin from 1 to 2 in the call to <b>adply</b> below.<br />
<br />
<pre>## box-and-whisker plot of FFT power spectrum by frequency/period
## vary these parameters -- r as per logistic map
## see http://en.wikipedia.org/wiki/Logistic_map for details
logmap = function(n, r, x0=0.5) {
ret = rep(x0,n)
for (ii in 2:n) {
ret[ii] = r*ret[ii-1]*(1-ret[ii-1])
}
return(ret)
}
## modify these for interesting behavior differences
rfrom = 3.4
rto = 3.7
rsteps=200
seqlen=1e4
require(plyr)
require(ggplot2)
mytimeseries = aaply(seq(from=rfrom, to=rto, length.out=rsteps), 1,
function(x) {
logmap(seqlen, x)
})
## you can plug in any array here for mytimeseries
## each row is a timeseries
## for series by column, change the margin from 1 to 2, below
logspec = adply( mytimeseries, 1, function(x) {
## change spec.pgram parameters as per goals
ret = spectrum(x, taper=0, demean=T, pad=0, fast=T, plot=F)
return( data.frame(freq=ret$freq, spec=ret$spec))
})
## boxplot stats at each unique frequency
logspec.bw = ddply(logspec, 'freq', function(x) {
ret = boxplot.stats(log10(x$spec));
## sometimes $out is empty, bind NA to keep dimensions correct
return(data.frame(t(10^ret$stats), out=c(NA,10^ret$out)))
})
## plot -- outliers are dots
## be careful with spacing of hard returns -- ggplot has odd rules
## as close to default "spectrum" plot as possible
ggplot(logspec.bw, aes(x=1/(freq))) + geom_ribbon(aes(ymin=X1,
ymax=X5), alpha=0.35, fill='green') +
geom_ribbon(aes(ymin=X2, ymax=X4), alpha=0.5, fill='blue') +
geom_line(aes(y=X3)) +
geom_point(aes(y=out), color='darkred') +
scale_x_continuous( trans=c('inverse'), name='Period') +
scale_y_continuous(name='Power', trans='log10')
</pre><br />
<br />
Following my nose further on the logistic map example, I also used the <b>animate</b> package to make a movie that walks through the power spectrum of the map for a range of r values. Maybe one day I'll set this to music and post it to Youtube! Though I think some strategic speeding up and slowing down in important parameter regions is warranted for a truly epic tour of the logistic map's power spectrum.<br />
<br />
<pre>require(animation)
saveVideo({
ani.options(interval = 1/5, ani.type='png', ani.width=1600, ani.height=1200)
for (ii in c(seq(from=3.2, to=3.5, length.out=200), seq(from=3.5, to=4, length.out=1200))) {
spectrum(logmap(1e4, ii), taper=0, demean=T, pad=0, fast=T, sub=round(ii, 3))
}
}, video.name = "logmap2.mp4", other.opts = "-b 600k")
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-45217636032646054912011-10-06T16:49:00.001-07:002011-10-06T16:49:58.702-07:00Bat CountryI've spent a lot of time thinking about and using R's <b>spectrum()</b> function and the <a href="http://en.wikipedia.org/wiki/Fast_Fourier_transform">Fast Fourier Transform</a> (FFT) in the last 5+ years. Lately, they've begun to remind me a little of a <a href="http://en.wikipedia.org/wiki/Theremin#Operating_principles">Theremin</a>: simple to use, difficult to master.<br />
<br />
While prepping a figure for the <a href="http://addictedtor.free.fr/graphiques/">R Graph Gallery</a> (which is awesome, by the way), I ran across this curious example -- its striking visual appearance was definitely <b>not</b> what I was expecting. <br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-tvmt2B3r29s/To4zzA8NmbI/AAAAAAAABCk/GZIIk1aGtCo/s1600/batcountry.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="400" width="400" src="http://2.bp.blogspot.com/-tvmt2B3r29s/To4zzA8NmbI/AAAAAAAABCk/GZIIk1aGtCo/s400/batcountry.png" /></a></div><br />
I decided to use the <a href="http://en.wikipedia.org/wiki/Logistic_map#Behavior_dependent_on_r">logistic map</a> to generate an ensemble of time series with a range of frequencies. The logistic map goes through characteristic period doubling, and then exhibits broad-spectrum "noise" at the onset of chaos. And it's <b>really</b> easy to compute. So, I can tune my time series to have a range of frequency distributions. <br />
<br />
In this example, I'm in the periodic regime (r=3.5, period 4) . So why am I getting this Batman motif?<br />
<br />
Actually, it's weirdly simple. This is an artifact of the (default) tapering that R's <b>spectrum()</b> function applies to the series before the FFT is computed. In theory, the Fourier Transform assumes an infinite-length sequence, while in practice the FFT assumes the series is circular -- in essence, gluing the beginning and end of the series together to make a loop. If the series is not perfectly periodic, though, this gluing introduces a sharp discontinuity. Tapering is typically applied to reduce or control this discontinuity, so that both ends gently decline to zero. When the series is genuinely periodic, though, this does some <b>weird</b> effects. In essence, the taper <a href="http://en.wikipedia.org/wiki/Transfer_function">transfers power</a> from the fundamental frequencies in a curious way. Note that the fundamental period is 4, as seen by the peak at frequency 1/4, with a strong harmonic at period 2, frequency 1/2.<br />
<br />
Zero-padding has a similar effect. If your series is zero (or relatively small) at both ends, then you could glue it together without introducing discontinuities, but peaks near the beginning and the end will still be seen as near by, resulting in spurious peaks. In any case, <b>R</b> pads your series by default (fast=T) unless the length of the series is highly composite (i.e. it can be factored into many small divisors). Below, it turns out that a series of length <b>1e5</b> is composite enough that <b>R</b> doesn't pad it. <b>1e5</b>? Totally different results.<br />
<br />
If I understand this correctly, both padding and tapering are attempts to reduce or control <a href="http://en.wikipedia.org/wiki/Spectral_leakage">spectral leakage</a>. We've constructed an example where no spectral leakage is expected from the original time series, which closely fits the assumptions of the FFT without modification.<br />
<br />
The moral? Know thy functions (and their default parameter values, and their effects)!<br />
<br />
<pre>## using R --
## logistic map function:
logmap = function(n, r, x0=0.5) {
ret = rep(x0,n)
for (ii in 2:n) {
ret[ii] = r*ret[ii-1]*(1-ret[ii-1])
}
return(ret)
}
## compare:
spectrum(logmap(1e5, 3.5), main='R tapers to 0.1 by default')
spectrum(logmap(1e5, 3.5), taper=0, main='No taper')
spectrum(logmap(1e5+1, 3.5), taper=0, main='Non-composite length, fast=T pads by default')
spectrum(logmap(1e5, 3.5), taper=0.5, sub='Lots of padding approximately recovers "correct" results')
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0tag:blogger.com,1999:blog-4695461192314112143.post-28605320053878761732011-06-18T20:49:00.014-07:002012-02-29T02:53:05.517-07:00Efficient loops in R -- the complexity versus speed trade-offI've <a href="http://helmingstay.blogspot.com/2011/05/problems-with-plyr-memorycomplexity.html">written before</a> about the up- and downsides of the <a href="http://cran.r-project.org/web/packages/plyr/index.html">plyr</a> package -- I love it's simplicity, but it can't be mindlessly applied, no pun intended. This week, I started building a agent-based model for a large population, and I figured I'd use something like a binomial per-timestep birth-death process for between-agent connections. <br />
<br />
My ballpark estimate was 1e6 agents using a weekly timestep for about 50 years. This is a stochastic model, so I'd like to replicate it at least 100 times. So, I'll need at least 50*50*100 = 250,000 steps. I figured I'd be happy if I could get my step runtimes down to ~1 second -- dividing the 100 runs over 4 cores, this would give a total runtime of ~17.5 hours. Not short, but not a problem.<br />
<br />
At first, I was disheartened to see the runtime of my prototype step function stretching into the minutes. What's going on? Well, I'd used <b>plyr</b> in a <b>very</b> inappropriate way -- for a very large loop. I began to investigate, and discovered that writing an aesthetically unappealing <b>while</b> gave me a 30+x speed-up. <br />
<br />
All of which got me thinking -- how expensive are loops and function calls in R? Next week I'm leading a tutorial in R and C++ using the wonderful <a href="http://cran.r-project.org/web/packages/Rcpp/index.html">Rcpp</a> and <a href="http://cran.r-project.org/web/packages/inline/index.html">inline</a> packages here at <a href="http://tuvalu.santafe.edu/events/workshops/index.php/Complex_Systems_Summer_School_2011">Santa Fe Institute's 2011 Complex Systems Summer School</a>. Might this make a nice example?<br />
<br />
It does, and in spades. Below are the test functions, and you can see that the code complexity increases somewhat from one to the other, but outcomes are identical, again with 30+x speedup for each subsequent case. Here, I'm using the native <a href="http://cran.r-project.org/doc/manuals/R-exts.html#The-R-API">R API</a>. I also tested using Rcpp to import <b>rbinom()</b>, but that ended up taking twice as long as the naive <b>while</b> loop.<br />
<br />
So, the moral of the story seems to be that if you can write a long loop in pure C++, it's a really easy win. <br />
<br />
<b>Note</b> -- The <b>as< double >(y);</b> in <b>src</b> below doesn't seem to copy-and-paste correctly for some reason. If <b>testfun2</b> doesn't compile, check to make sure this bit pasted correctly. <br />
<br />
<br />
<h3>The pure R function definitions</h3><pre>## use aaply -- the simplest code
require(plyr)
testfun0 <- function(x, y) aaply(x, 1, function(x) rbinom(1, x, y)
## rewrite the above as an explicit loop
testfun1 = function(nrep, x, y) {
## loop index
ii<-0;
## result vector
ret<-rep(0, nrep);
while(ii < nrep) {
ii<-ii+1;
## how many successes for each element of bb?
ret[ii] <- rbinom(1, x[ii], y)
}
return(ret)
}
</pre>
<h3>Rcpp function definitions (with lots of comments)</h3><pre>## define source code as string, pass to inline
src <- '
// Cast input parameters to correct types
// clone prevents direct modification of the input object
IntegerVector tmp(clone(x));
double rate = as< double >(y);
// IntegerVector inherits from STL vector, so we get standard STL methods
int tmpsize = tmp.size();
// initialize RNG from R set.seed() call
RNGScope scope;
// loop
for (int ii =0; ii < tmpsize; ii++) {
// Rf_rbinom is in the R api
// For more details, See Rcpp-quickref.pdf, "Random Functions"
// also, Rcpp provides "safe" accessors via, e.g., tmp(ii)
tmp(ii) = Rf_rbinom(tmp(ii), rate);
};
// tmp points to a native R object, so we can return it as-is
// if we wanted to return, e.g., the double rate, use:
// return wrap(rate)
return tmp;
'
require(inline)
## compile the function, inspect the process with verbose=T
testfun2 = cxxfunction(signature(x='integer', y='numeric'), src, plugin='Rcpp', verbose=T)
## timing function
ps <- function(x) print(system.time(x))
</pre>
<h3>The tests</h3><pre>## Input vector
bb <- rbinom(1e5, 20, 0.5)
## test each case for 2 different loop lengths
for( nrep in c(1e4, 1e5)){
## initialize RNG each time with the same seed
## plyr
set.seed(20); ps(cc0<- testfun0(nrep[1:nrep], bb, 0.1))
## explicit while loop
set.seed(20); ps(cc1<- testfun1(nrep, bb, 0.1))
## Rcpp
set.seed(20); ps(cc2<- testfun2(bb[1:nrep], 0.1))
print(all.equal(cc1, cc2))
}
## output
user system elapsed
3.392 0.020 3.417
user system elapsed
0.116 0.000 0.119
user system elapsed
0.000 0.000 0.002
[1] TRUE
user system elapsed
37.534 0.064 37.601
user system elapsed
1.228 0.000 1.230
user system elapsed
0.020 0.000 0.021
[1] TRUE
</pre>
<h3>Postlude</h3>After posting this to the very friendly <a href="http://lists.r-forge.r-project.org/pipermail/rcpp-devel/">Rcpp-devel</a> mailing list, I got an in-depth reply from Douglas Bates pointing out that, in this case, the performance of vanilla <b>apply()</b> beats the <b>while</b> loop by a narrow margin. He also gives an interesting example of how to use an STL template and an iterator to achieve the same result. I admit that templates are still near-magic to me, and for now I prefer the clarity of the above. Still, if you're curious, this should whet your appetite for some of the wonders of Rcpp.
<pre>## Using Rcpp, inline and STL algorithms and containers
## use the std::binary_function templated struct
inc <- '
struct Rb : std::binary_function< double, double, double > {
double operator() (double size, double prob) const { return ::Rf_rbinom(size, prob); }
};
'
## define source code as a character vector, pass to inline
src <- '
NumericVector sz(x);
RNGScope scope;
return NumericVector::import_transform(sz.begin(), sz.end(), std::bind2nd(Rb(), as< double >(y)));
'
## compile the function, inspect the process with verbose=TRUE
f5 <- cxxfunction(signature(x='numeric', y='numeric'), src, 'Rcpp', inc, verbose=TRUE)
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com10tag:blogger.com,1999:blog-4695461192314112143.post-40124924195054899912011-06-11T19:29:00.001-07:002012-02-29T02:54:25.108-07:00The importance of being unoriginal (and befriending google)<h3>In search of bin counts</h3>I look at histograms and density functions of my data in R on a regular basis. I have some idea of the algorithms behind these, but I've never had any reason to go under the hood until now. Lately, I've been looking using the bin counts for things like Shannon entropy ( in the very nice <a href="http://cran.r-project.org/web/packages/entropy">entropy</a> package. I figured that binning and counting data would either be supported via a native, dedicated R package, or quite simple to code. Not finding the former ( <texttt>base graphics hist()</texttt> uses <texttt>.Call("bincounts")</texttt>, which appears undocumented and has a boatload of arguments ), I naively failed to search for a package and coded up the following.<br />
<br />
<pre>myhist = function(x, dig=3) {
x=trunc(x, digits=dig);
## x=round(x, digits=dig);
aa = bb = seq(0,1,1/10^dig);
for (ii in 1:length(aa)) {
aa[ii] = sum(x==aa[ii])
};
return(cbind(bin=bb, dens=aa/length(x)))
}
## random variates
test = sort(runif(1e4))
get1 = myhist(test)
</pre><br />
<h3>Trouble in paradise</h3>Truncate the data to a specified precision, and count how many are in each bin. Well, first I tried <texttt>round(x)</texttt> instead of <texttt>trunc(x)</texttt>, which sorta makes sense but gives results that I still don't understand. On the other hand, <texttt>trunc(x)</texttt> doesn't take a digits argument? WTF? Of course, I could use <texttt>sprintf(x)</texttt> to make a character of known precision and convert back to numeric, but string-handling is waaaaaay too much computational overhead. Like towing a kid's red wagon with a landrover...<br />
<br />
<h3>Dear Google...</h3>An hour of irritation and confusion later, I <a href="http://www.google.com/search?q=cran+bin+counts">ask google</a> and, small wonder, the second search result links to the <a href="http://cran.r-project.org/web/packages/ash">ash</a> package that contains said tool. And it runs somewhere between 100 and 1,000 times faster. It doesn't return the bin boundaries by default, but it's good enough for a quick-and-dirty empirical probability mass distribution.<br />
<br />
To be fair, there's something to be said for cooking up a simple solution to a simple problem, and then realizing that, for one reason or another, the problem isn't quite as simple as one first thought. On the other hand, sometimes we just want answers. When that's the case, asking google is a pretty good bet. <br />
<br />
<pre>## their method
require(ash)
get2 = bin1(test, c(0,1), 1e3+1)$nc
</pre>xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com3tag:blogger.com,1999:blog-4695461192314112143.post-24581173968789085952011-05-10T03:42:00.001-07:002012-02-29T02:51:31.713-07:00Problems with plyr -- the memory/complexity trade-off<h3>Two types of R users</h3><br />
My overwhelming impression from UseR 2010 is that, generally speaking, there are 2 types of regular R users -- those who have heard and are made uncomfortable by the idea of the <b>*apply()</b> functions, and those who really <em>get</em> it. In UNM R programming group that I've been leading for about a year now, I've really tried to get people over the hump and into the second group. Once there, folks seem to really appreciate the amazing power of vectorization in R, and begin to enjoy writing code. The conceptual clarity of:<br />
<br />
<pre>mymatrix = matrix(1:10, nrow=10)
apply(mymatrix, 1, sum)
apply(mymatrix, 2, sum)
</pre><br />
over the clunky:<br />
<br />
<pre>rowSums(mymatrix)
colSums(mymatrix)
</pre><br />
may not be immediately apparent. Eventually, though, folks go searching for things like <b>rowMedian()</b> and become frustrated that R doesn't "have all the functions" that they need. Well, R does, once you grok <b>apply()</b>.<br />
<br />
<h3>Hadley's Magic Brainchild </h3>In the last year, I've had some serious <b>ahah!</b> moments with <a href="http://had.co.nz/plyr/">plyr</a>. Lately, I've added the <a href="http://had.co.nz/reshape/">reshape</a> package to the mix to achieve some serious R 1-liner Zen. Multidimensional arrays to long-form dataframes?<br />
<br />
<pre>myarray = array(0, dim=c(3,4,5,6,10), dimnames=list(a=1:3, b=1:4, c=1:5, d=letters[1:6], e=LETTERS[1:10]))
melt( adply(myarray, 1:4, I))
</pre><br />
Sure. No problem. It's really that easy?!<br />
<p><br />
An unexpeceted bonus? This way of thinking lends itself nicely to thinking of explicit parallelism. If you can describe the problem as a single recipe that's done for each "piece" of a whole, then you're one step away from using all the cores on your machine to solve a problem with <a href="http://cran.r-project.org/web/packages/foreach/index.html">foreach</a>. Need to apply some sort of long-running analysis to each element of a list, and want to return the results as a list? Why write a <b>for</b> loop when you can do:<br />
<br />
<pre>## prep
## warning -- this may not be fast
## use smaller matrices for smaller machines!
nn = 2^20
rr = 2^10
mylist = list(a=matrix(rnorm(nn), nrow=rr), b=matrix(rnorm(nn), nrow=rr), c=matrix(rnorm(nn), nrow=rr))
## analysis
myeigs = foreach( ii = iter(mylist), .combine=c) %do% { print('ping'); return(list(eigen(ii)))}
</pre><br />
and then, once it works, change <b>%do%</b> to <b>%dopar%</b>, add the following, and you're off to the races!<br />
<br />
<pre>require(multicore)
require(doMC)
## use the appropriate # of cores, of course
registerDoMC(cores=4)
myeigs1 = foreach( ii = iter(mylist), .combine=c) %dopar% { print('ping'); return(list(eigen(ii)))}
</pre><br />
Compared to something like <b>llply</b>, your <b>dimnames</b> don't automatically propagate to the results, but I think this is still pretty amazing. Complete control <em>and</em> simple parallelism.<br />
<br />
Debugging with <b>%dopar%</b> is tricky, of course, because there are separate stacks for each item (i think), and messages (such as calls to <b>print()</b>) don't return to the console as you might expect them to. So, when in doubt, back off to <b>%do%</b>.<br />
<br />
<h3>What could possibly go wrong?</h3>The only problem with all of this is, when tasks are embarassingly parallel, that data also becomes embarassingly parallel to point where it no longer fits into memory. Thus, I returned today to a long-running bootstrap computation to find R consuming ~7GB RAM, 5+ GB swap, and this message:<br />
<br />
<pre>Error: cannot allocate vector of size 818.6 Mb
Enter a frame number, or 0 to exit
...
4: mk.bootrep(zz, 75)
5: mk.boot.R#165: aaply(tmpreps, .margins = 1, function(x) {
6: laply(.data = pieces, .fun = .fun, ..., .progress = .progress, .drop = .dro
7: list_to_array(res, attr(.data, "split_labels"), .drop)
8: unname(unlist(res))
</pre><br />
What's happening is that <b>plyr</b> is trying to do everything at once. As anyone who's used <b>grep</b> can tell you, doing one row at a time, or <em>streaming</em> data is often a much better idea. I got the intended results from above by pre-allocating an array and writing each item of my results list into the array in a loop in seconds, and barely broke 3 GB of RAM usage.<br />
<br />
Now, nothing here is really news. The dangers of "Growing Objects" is covered in Circle 2 of Burns Statistics' wonderful <a href="http://www.burns-stat.com/pages/Tutor/R_inferno.pdf">R Inferno</a>. Still, <b>plyr</b> strikes me as an interesting case where reducing conceptual complexity can lead to a rather steep increase in computational complexity. And the most interesting thing of all is that it happens quite suddenly above a certain threshold.<br />
<br />
<h3>Parting thoughts</h3>I wonder if there's any major barrier to a <b>stream=TRUE</b> argument to the <b>plyr</b> functions -- I haven't thought about it too much, but imagine that you'd also need a finalizer function to prepare the return object to be written/streamed into. At what point is it easier to just do by hand with a <b>for</b> loop? <br />
<br />
Honestly, I don't know the answer. I don't do too many things that break <b>plyr</b>, but I've learned how important it is to understand when I'm likely exceed its limits. <b>.variables</b> in <b>ddply</b> is another one that I've learned to be careful with. If, after subdividing your input data.frame, you end up with 1e5 or 1e6 pieces, things start to break down pretty fast.<br />
<br />
Honestly, I love writing solutions with the <b>*apply()</b> family and the <b>ddply</b> functions. I think it makes code cleaner and more logical. Watching the light of this dawn in other R users' eyes is truly exciting. Yet it pays to remember that, like all things, it has its limits. In this case, one seems to reach the limits suddenly and harshly.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com1tag:blogger.com,1999:blog-4695461192314112143.post-24233675029194951372011-03-29T04:16:00.002-07:002012-02-29T02:55:31.345-07:00Why PKI matters.<a href="https://www.eff.org/deeplinks/2011/03/iranian-hackers-obtain-fraudulent-https">This article</a> landed in my inbox this week in a newsletter from the EFF. I usually don't read them, but the term "meltdown" caught my eye, what with all the nuke new this month. They also managed to work in "too big to fail", and neither reference was hyperbolic. The internet depends on a level of trust, and (surprise) there are people working to co-opt that trust.<br />
<br />
A number of my friends think I'm a little weird for using <emph>real</emph> passwords, for not sharing them, etc. But I read Cryptonomicon and friends; I have at least one box exposed to the real, honest-to-god, jungle-out-there internet; I have a healthy fear of all the things that could go wrong. As part of my <a href="http://math.unm.edu/~hwearing/researchgroup/index.html">lab group's</a> data entry project, I recently registered my first SSL certificate from <a href="http://www.startssl.com/">http://www.startssl.com/</a> (free!), and learned quite a bit about PKI in the process.<br />
<br />
Still, reading an article like this really drives home both the complexities and importance of the global PKI system. Trust is difficult when it's strung across the globe on fiber optics cables, and enforced by our inability to quickly factor very large numbers. But old-school techniques of impersonation and breaking and entering will always be with us. I may trust google.com, but how do I know that it's actually them? PKI.<br />
<br />
How cool is that?<br />
Very.xianhttp://www.blogger.com/profile/01991290472959345882noreply@blogger.com0