Sunday, November 24, 2013

Just for fun: attractors in R

I have a borderline unhealthy obsession with attractors. I thought I got it out of my system, but here we are. For whatever reason, I felt like making some in R.

You can find the R code here. It uses the attractor function to define density in a matrix, which is how often a given point gets hit. Then we plot the log scaled version of that using image().

Be warned, with a lot of points it can be quite slow. If you play around you might want to drop the n value to 10,000 or so.

Some pictures

Dejong


Clifford



If you like this sort of thing, I made a realtime interactive dejong attractor, using a kinect, openFrameworks and a lot of GLSL shaders. You can see a clip of it below, the source is also on my github if you are so inclined.


Attractor from dizzy pete on Vimeo.

You can find me on twitter and G+

10 comments:

  1. Nice. Without modifications you can speed up the processing time ~2 times by compiling the functions:

    library("compiler")
    setCompilerOptions(suppressAll=TRUE)
    map <- cmpfun(map)
    dejong <- cmpfun(dejong)

    ReplyDelete
  2. This comment has been removed by the author.

    ReplyDelete
  3. This comment has been removed by the author.

    ReplyDelete
  4. With vectorization you can speed it up 20-40 times:

    mapxy <- function(x, y, xmin, xmax, ymin=xmin, ymax=xmax) {
    row <- round( (x - xmin) / (xmax - xmin) * (width - 1) + 1 )
    col <- round( (y - ymin) / (ymax - ymin) * (height - 1) + 1 )
    (col-1)*height + row
    }

    dejong2 <- function(x, y) {
    for (i in 1:npoints) {
    xn <- sin(a * y) - cos(b * x)
    yn <- sin(c * x) - cos(d * y)
    idxs <- mapxy(xn, yn, -2, 2)
    mat[] <<- mat[] + tabulate(idxs, nbins=width*height)
    x <- xn
    y <- yn
    }
    }

    clifford2 <- function(x, y) {
    for (i in 1:npoints) {
    xn <- sin(a * y) + c * cos(a * x)
    yn <- sin(b * x) + d * cos(b * y)
    idxs <- mapxy(xn, yn, -abs(c)-1, abs(c)+1, -abs(d)-1, abs(d)+1)
    mat[] <<- mat[] + tabulate(idxs, nbins=width*height)
    x <- xn
    y <- yn
    }
    }

    setCompilerOptions(suppressAll=TRUE)
    mapxy <- cmpfun(mapxy)
    dejong2 <- cmpfun(dejong2)
    clifford2 <- cmpfun(clifford2)

    which you call as:

    dejong2(rsamp[,1], rsamp[,2])
    clifford2(rsamp[,1], rsamp[,2])

    ReplyDelete
    Replies
    1. hey thanks, that is really cool. It didn't even occur to me to do that. I will add it to the script

      Delete
    2. Yet faster is to do:

      mapxy <- function(x, y, xmin, xmax, ymin=xmin, ymax=xmax) {
      sx <- (width - 1) / (xmax - xmin)
      sy <- (height - 1) / (ymax - ymin)
      row0 <- round( sx * (x - xmin) )
      col0 <- round( sy * (y - ymin) )
      col0 * height + row0 + 1
      }

      dejong3 <- function(x, y) {
      nidxs <- length(mat)
      counts <- integer(length=nidxs)
      for (i in 1:npoints) {
      xt <- sin(a * y) - cos(b * x)
      y <- sin(c * x) - cos(d * y)
      x <- xt
      idxs <- mapxy(x, y, -2, 2)
      counts <- counts + tabulate(idxs, nbins=nidxs)
      }
      mat <<- mat + counts
      }

      clifford3 <- function(x, y) {
      ac <- abs(c)+1
      ad <- abs(d)+1
      nidxs <- length(mat)
      counts <- integer(length=nidxs)
      for (i in 1:npoints) {
      xt <- sin(a * y) + c * cos(a * x)
      y <- sin(b * x) + d * cos(b * y)
      x <- xt
      idxs <- mapxy(x, y, -ac, ac, -ad, ad)
      counts <- counts + tabulate(idxs, nbins=nidxs)
      }
      mat <<- mat + counts
      }

      Delete
    3. Hey thanks Henrik, have updated the script on github. It really is amazing how much faster it is. Apologies that it took so long to update, just caught up with some other things.

      Thanks again!

      Delete
  5. @ Henrik: Extremely cool use of matrix addition using mat[ ]! Wondering if Rcpp will be even faster, but I don't think much...

    ReplyDelete
  6. @ pete: One can also play around using erosion of small value points, .i.e.

    QUANT <- quantile(mat, 0.5)
    mat[mat <= QUANT] <- 0
    dens <- log(mat + 1)/round(log(max(mat)))

    ReplyDelete