Showing posts with label R. Show all posts
Showing posts with label R. Show all posts

Thursday, August 18, 2011

Do older SOers use fewer words?

On StackOverflow, to posters with more experience ask their questions in fewer words?

No. There's no visible difference:

Chars of non-code:

so.rst-weaver/description.png

Chars of code:

so.rst-weaver/code.png

The data comes from the super-handy StackOverflow API, which was retrieved using wget and then parsed using rjson and XML.

First read in and parse the JSON:

     so.R 

  1  library(rjson)
  2  library(XML)
  3  library(ggplot2)
  4  library(plyr)
  5  
  6  read.qs = function(path) {
  7      fromJSON(file = path)$questions
  8  }
  9  
 10  questions = do.call(c,
 11      lapply(c('page-1.json', 'page-2.json', 'page-3.json'),
 12          read.qs
 13      )
 14  )

Then for each one parse the HTML and look for <pre> and <p> tags:

     so.R (cont)

 15  Table = ldply(questions, function(q) {
 16      body.text = sprintf('<body>%s</body>', q$body)
 17      body = htmlParse(body.text)
 18  
 19      description = tot.length.of(body, '//p//text()')
 20      code = tot.length.of(body, '//pre//text()')
 21  
 22      rep = q$owner$reputation
 23  
 24      data.frame(
 25          rep, description, code
 26      )
 27  })

(where tot.length.of is:

     so.R (cont)

 28  tot.length.of = function(doc, query) {
 29      parts = xpathApply(doc, query, xmlValue)
 30      text = paste(parts, collapse='')
 31      nchar(text)
 32  }

)

Then make the plots:

     so.R (cont)

 33  png('description.png')
 34  print(ggplot(data=Table)
 35      + geom_point(aes(rep, description))
 36      + scale_x_log10()
 37      + scale_y_log10()
 38      + xlab('Rep')
 39      + ylab('Verbosity')
 40  )
 41  dev.off()
 42  
 43  png('code.png')
 44  print(ggplot(data=Table)
 45      + geom_point(aes(rep, code))
 46      + scale_x_log10()
 47      + scale_y_log10()
 48      + xlab('Rep')
 49      + ylab('Verbosity')
 50  )
 51  dev.off()

$ Rscript so.R >/dev/null 2>&1

Monday, July 11, 2011

You can scrap it and write something better but let me keep R ;)

Ross Ikaha (via Xi'an -- thanks ;) ) gives a nice example to show why R is basically impossible to optimize:

> f = function() {

> if (runif(1) > 0.5) {

> x = 10

> }

> x

> }

x in the last expression could be a local or a global, and this won't be known until runtime!

That's pretty good.

I wonder if this will always return 10:

> f = function() {

> a = 10

> a

> }

> f()

[1] 10

Maybe we could optimize it out? No:

> regular.equals <- `=`

> `=` = function(...) {

> with(parent.frame(), a <- 11)

> }

> f()

[1] 11

> `=` <- regular.equals

OK but assigning to `=` is really cheating. This is cheating too:

> regular.brace = `{`

> `{` = function(...) 11

> f()

[1] 11

> `{` = regular.brace

If `=` and `{` weren't available for assignment, is there any way to make f not return 10? Well there's this, but it's hardly "making f return 10":

> body(f) = quote(11)

> f()

[1] 11

I can't think of any other way to attack it; so if you block all those I guess we can declare this function safe ;)

How about this one:

> g = function(x, y) {

> x + y

> }

Could we safely say that while g() is being called, f() will never be called at the same time (ie so that we could overlay their stacks?)

No, not even close. There are a million ways to break that:

> g(f(), f())

[1] 22

(Lazy evaluation -- f() is called while executing g())

> normal.plus = `+`

> `+` = function(a, b) UseMethod("+")

> `+.numeric` = function(a, b) { f() * a * b }

> g(1, 2)

[1] 22

> `+` = normal.plus

(Generic `+` is not uncommon).

Well at least can we say g() will always evaluate both arguments? No:

> g(return(2), print("hi"))

[1] 2

Which can make code that looks like an error be valid:

> g(return(2), not.defined)

[1] 2

If you assume g() is strict you will knock out some perfectly legitimate R code :)

Likewise an explicit call to return() doesn't have to actually do anything:

> g = function(x) {

> 3

> }

> g(return(2))

[1] 3

I also find these hacks delightful:

> delayedAssign('b', {b <- 7; 8})

> a <- b

> a

[1] 8

> b

[1] 7

> x = '@'

> up = function() {

> x. <- x

> x <<- intToUtf8(utf8ToInt(x)+1L)

> delayedAssign(x, up(), ass=.GlobalEnv)

> x.

> }

> up()

[1] "@"

> A

[1] "A"

> B

[1] "B"

> C

[1] "C"

> D

[1] "D"

> backw = function(...) {

> args = as.list(substitute(list(...)))[-1L]

> env = parent.frame()

>

> for (arg in rev(args)) {

> res <- eval(arg, env)

> }

> res

> }

> environment(backw) = baseenv()

> `{` = backw

> fib = function(n) {

> a

> for (i in 1:n) {

> a = b - a

> b = a + b

> }

> b = 1

> a = 1

> }

> fib(8)

[1] 34

> lex.scope = function(...) {

> args = as.list(substitute(list(...)))[-1L]

> parent = parent.frame()

> env = new.env(parent=parent)

>

> for (arg in args) {

> res <- eval(arg, env)

> }

> res

> }

> environment(lex.scope) = baseenv()

> `{` = lex.scope

> {

> x = 2

> {

> x = 3

> }

> x

> }

[1] 2

(Hey, that one could actually be useful).

> evanescent = function(x, v) {

> force(v)

> env = parent.frame()

> name = as.character(substitute(x))

>

> delayedAssign(name, {

> env[[name]] = alist(x=)$x

> v

> }, ass = env

> )

> }

> environment(evanescent) = baseenv()

> `=` = evanescent

> a = 2

> a

[1] 2

> a

Error in eval(expr, envir, enclos) :

argument "a" is missing, with no default

My thoughts on all this: R is a wonderful and flexible language, but it is entirely not the right language for intense numerical calculation: it is simply too flexible to be fast. And unfortunately whenever you're trying some algorithm that hasn't been done before, it will necessarily not already be written in C, and you find yourself having to make a choice: at what point is it too slow to keep in R, and the best thing to do is to give up and rewrite it in C?

The thing is, as an interactive environment, R's flexibility really gains something (think how you'd implement with() in Python). I personally am OK with the current system of using R for high level code and ducking to C for tight loops; it's better than trying to make one language be good at 2 different things.

Sunday, July 3, 2011

Reverse Iteration

Time to horrify some people.

First let's include the code we wrote last time,

> source("pretend.R")

and the dependency-tracking environment it creates will be used to run all the following examples.

Let's look at, I don't know, I'm just trying to demonstrate a language feature so uh... band-pass filtering Gaussian noise.

Here's some noise:

> t = seq(0, 10, len=40)

> x = rnorm(length(t))

> plot(t, x, type='l')

And filter it:

> X = rfft(x)

> f = 1/length(t)/(t[2]-t[1]) * (seq_along(X)-1)

> s = 2*pi*1i*f

> B = .2

> f0 = .5

> H = B/(s**2 + B*s + (2*pi*f0)**2)

> Y = H*X

> y = irfft(Y)

> min(f)

[1] 0

> max(f)

[1] 1.95

So that it looks like:

> plot(t, y, type='l')

We've made it look wavier. We could decrease the bandwidth and make it cleaner:

> B = .1

> plot(t, y, type='l')

We could unrwap the phase and find the (what we already know to be) the frequency of the wave we just made:

> y.a = analytic(y)

> phase = unwrap(Arg(y.a))

> f0.est = qr.solve(

> cbind(t, 1),

> cbind(phase)

> )[[1]] / (2*pi)

> f0.est

[1] 0.4818572

We can try various bandwidths and see what estimates we get for the center:

> Bs = seq(0.1, 0.3, len=12)

> f0.ests = sapply(Bs, function(B.) {

> B <<- B.

> f0.est

> })

> f0.ests

[1] 0.4818572 0.4814956 0.4811829 0.4809197 0.4807073 0.4805482 0.4804464
[8] 0.4804080 0.4804399 0.4805393 0.4806539 0.3708171

> plot(Bs, f0.ests)

And make the same plot for a different center:

> f0 = 1.0

> f0.ests

[1] 0.9406669 0.9400220 0.9393042 0.9385253 0.9376873 0.9367848 0.9358192
[8] 0.9348181 0.9338298 0.9328914 0.9320139 0.8671673

> plot(Bs, f0.ests)

Now the interesting thing about all that we've done here is that nowhere is the interesting code wrapped in either functions or loops: it's not wrapped in anything at all. There are loops, but they occurr after. Code such as

> y.a = analytic(y)

> phase = unwrap(Arg(y.a))

> f0.est = qr.solve(cbind(t, 1), cbind(phase))[[1]]/(2 * pi)

was written entirely without repetition in mind.

I have to admit that I did cheat a little bit. If you go back and look at the code you'll notice that each variable is assigned exactly once. Hence the hiding of logic like

> analytic = function(u) {

> n = length(u)

> m = n/2 + 1

> U = fft(u)/length(u)

> U[(m + 1):n] = 0

> fft(U, inv = T)

> }

into functions, because this is awkward to write without assigning U twice. Because in the main program each variable is only assigned once, when a modification is made, it is trivial to decide which other variables must be updated and how.

Sunday, June 26, 2011

RObjectTables are AWESOME

Why isn't everyone using the RObjectTables package? This is the best thing ever!

Here's the basic idea of RObjectTables: An environment is an object where you can lookup names and associate them with values. And in particular its where you look up variables used in an expression. But there's no reason you can't take any other object that associates names with values (data.frame, list, SQL database, CSV file, filesystem, ...) and use that to lookup your variables.

R already has sortof this because you can use eval() and friends on lists and data.frames. But it's not extensible. The purpose of RObjectTables is to make it extensible.

In its current version RObjectTables is somewhat limited, because you can only attach() the created environments, not pass them to with() and such. But with extremely minor modifications this becomes possible. You can find a version with these changes on my GitHub site.

Now let me demonstrate how unbelievably useful this is. So naturally I will start with a useless example: an environment of only strings!

> library(RObjectTables)

> db = newRClosureTable(list(

> assign = function(name, value) {

> # Not used

> },

>

> get = function(name) {

> name

> },

>

> exists = function(name) {

> T

> },

>

> remove = function(name) {

> # Not used

> },

>

> objects = function() {

> # Not used

> }

> ))

> with(db, x)

[1] "x"

Muahahaha.

This environment, as you can see, is not terribly useful:

> try(

> with(db, x + y),

> silent=T

> )

> geterrmessage()

[1] "Error in x + y : non-numeric argument to binary operator\n"

Yes... using "+" (character) as a function... not so good. Also we don't have to feel so bad about not implementing assign() because with `<-` turning into "<-" there's not much we could assign.

But now let's do something really useful.

> reality = function() {

> parent = parent.frame()

> outer = new.env(parent = parent)

>

> formulas = list()

> valid = list()

> values = list()

>

> make.ref = function(name) {

> ref = list(

> db = self,

> name = name

> )

> class(ref) = 'ref'

> ref

> }

>

> dep.table = list()

> rdep.table = list()

>

> propagate = function(name) {

> if (is.null(valid[[name]]) || !valid[[name]]) {

> }

> else {

> for (ref in rdep.table[[name]]) {

> ref$db$reset(ref$name)

> }

> }

> }

>

> reset = function(name) {

> propagate(name)

> valid[[name]] <<- F

> }

>

> add.dep = function(name, ref) {

> dep.table[[name]] <<- c(list(ref), dep.table[[name]])

> }

>

> add.rdep = function(name, ref) {

> rdep.table[[name]] <<- c(list(ref), rdep.table[[name]])

> }

>

> del.rdep = function(name, ref) {

> rdep.table[[name]] <<- setdiff(rdep.table[[name]], list(ref))

> }

>

> ptr = newRClosureTable(list(

> assign = function(name, value) {

> force(value)

> formulas[[name]] <<- function() value

> reset(name)

> },

> get = function(name) {

> if (is.null(formulas[[name]])) {

> return(tryCatch(

> get(name, outer),

> error = function(e) getUnbound()

> ))

> }

>

> this.ref = make.ref(name)

> for (ref in working.refs) {

> ref$db$add.dep(ref$name, this.ref)

> }

>

> if (is.null(valid[[name]]) || !valid[[name]]) {

> old.deps <<- dep.table[[name]]

> dep.table[[name]] <<- list()

>

> working.refs <<- c(list(this.ref), working.refs)

>

> values[[name]] <<- formulas[[name]]()

>

> working.refs <<- working.refs[-1L]

>

> lost.deps = setdiff(old.deps, dep.table[[name]])

> for (ref in lost.deps) {

> ref$db$del.rdep(ref$name, this.ref)

> }

>

> gained.deps = setdiff(dep.table[[name]], old.deps)

> for (ref in gained.deps) {

> ref$db$add.rdep(ref$name, this.ref)

> }

>

> valid[[name]] <<- T

> }

>

> values[[name]]

> },

> exists = function(name) {

> !is.null(valid[[name]]) || exists(name, parent)

> },

> remove = function(name) {

> # TODO: this

> },

> objects = function(name) {

> names(valid)

> }

> ))

> class(ptr) = c('reality', class(ptr))

>

> attr(ptr, 'delayedAssign') = (

> function(name, promise) {

> formulas[[name]] <<- promise

> reset(name)

> promise

> }

> )

>

> self = list(

> ptr = ptr,

> reset = reset,

> add.dep = add.dep,

> add.rdep = add.rdep,

> del.rdep = del.rdep

> )

> class(self) = 'realptr'

>

> ptr

> }

> working.refs = list()

Why did I call it a "reality"? Well... uh... all the good names are taken. So think of it as storing the current version of... reality. What it does is keep track of a set of variables that depend on each other. Oh yeah and we need some way to assign into it:

> `$.reality` = function(r, name) {

> get(name, envir=r)

> }

> `$<-.reality` = function(r, name, value) {

> expr = substitute(value)

> env = parent.frame()

>

> attr(r, 'delayedAssign')(

> name, function() eval(expr, env)

> )

>

> r

> }

Now we can write:

> r = reality()

> r$x = 10

> r$y = sqrt(r$x)

> r$y

[1] 3.162278

> r$x = 20

> r$y

[1] 4.472136

It is the feature I have always been wishing for. Now if only we could get .GlobalEnv to do this... You might think of attach()ing it: don't do that. Trust me. I can't show you because it will crash Sweave and not produce any output, but that doesn't work yet (it would be nice if it did).

Luckily R has a useful feature that was never intended for this purpose (I find it is this way with most of R's finest features): the debugging browser. So let us...

> with.db = function(env, func1, func2) {

> if ('reality' %in% class(env)) {

> func1(env)

> }

> else {

> func2()

> }

> }

> reality.assign = function(name, rvalue) {

> env <- parent.frame()

> name <- substitute(name)

> srvalue <- substitute(rvalue)

>

> with.db(env,

> function(db) {

> name <- as.character(deparse(name))

> promise <- function() eval(srvalue, env)

>

> attr(db, 'delayedAssign')(name, promise)

> alist(x=)$x

> },

> function() {

> do.call(`<-`, list(name, rvalue), envir=env)

> }

> )

> }

> with(r, `=` <- reality.assign)

> run.interpreter = function(base.env) {

> options(browserNLdisabled = T)

> with(base.env, browser())

> }

Now I don't know how to feed an interactive debugging session into Sweave, so here is copy-pasted this feature in use:

Browse[1]> y = sqrt(x)

Browse[1]> for (i in 1:10) {
+ x = i
+ print(y)
+ }
[1] 1
[1] 1.414214
[1] 1.732051
[1] 2
[1] 2.236068
[1] 2.449490
[1] 2.645751
[1] 2.828427
[1] 3
[1] 3.162278
Browse[1]> b = a

Browse[1]> c = b

Browse[1]> a = 1

Browse[1]> c
> run.interpreter(r)
Called from: eval(expr, envir, enclos)
Browse[1]> y = sqrt(x)

Browse[1]> for (i in 1:10) {
+ x = i
+ print(y)
+ }
[1] 1
[1] 1.414214
[1] 1.732051
[1] 2
[1] 2.236068
[1] 2.449490
[1] 2.645751
[1] 2.828427
[1] 3
[1] 3.162278
Browse[1]> z = y**2

Browse[1]> z
[1] 10
Browse[1]> x = 1:3

Browse[1]> z
[1] 1 2 3
Browse[1]> b = a

Browse[1]> c = b

Browse[1]> a = 1

Browse[1]> c
> 

Haha "c" closed the debugger ;)

So far this actually seems stable enough for general use. I wouldn't trust it with data you care about, but I intend to make great use of it.