Plotting strange attractors using `datashader` with `reticulate`

strange_attractor.R

# Original blog post
# https://fronkonstin.com/2019/01/10/rcpp-camaron-de-la-isla-and-the-beauty-of-maths/

# For this example to work, you must have a conda environment with the datashader library
# See http://datashader.org/getting_started/index.html#installation for installation instructions

library(Rcpp)
library(reticulate)
library(dplyr)
use_condaenv("datashader", required = TRUE)

ds <- import("datashader")
pd <- import("pandas")



cppFunction('
  DataFrame createTrajectory(int n, double x0, double y0, 
  double a1, double a2, double a3, double a4, double a5, 
  double a6, double a7, double a8, double a9, double a10, 
  double a11, double a12, double a13, double a14) {
  // create the columns
  NumericVector x(n);
  NumericVector y(n);
  x[0]=x0;
  y[0]=y0;
  for(int i = 1; i < n; ++i) {
    x[i] = a1+a2*x[i-1]+ a3*y[i-1]+ a4*pow(fabs(x[i-1]), a5)+ a6*pow(fabs(y[i-1]), a7);
    y[i] = a8+a9*x[i-1]+ a10*y[i-1]+ a11*pow(fabs(x[i-1]), a12)+ a13*pow(fabs(y[i-1]), a14);
  }
  // return a new data frame
  return DataFrame::create(_["x"]= x, _["y"]= y);
}
')

{
  a1 <- -0.8
  a2 <- 0.4
  a3 <- -1.1
  a4 <- 0.5
  a5 <- -0.6
  a6 <- -0.1
  a7 <- -0.5
  a8 <- 0.8
  a9 <- 1.0
  a10 <- -0.3
  a11 <- -0.6
  a12 <- -0.3
  a13 <- -1.2
  a14 <- -0.3
}

library(scales)
# scales::trans
# scales::col_numeric()

par(mar = rep(0, 4), mai = rep(0, 4))
system.time({
  
  df <- 
    createTrajectory(1e7, 1, 1, a1, a2, a3, a4, a5, a6, 
                     a7, a8, a9, a10, a11, a12, a13, a14) %>% 
    filter(
      x > quantile(x, probs = 0.05, na.rm = TRUE), 
      x < quantile(x, probs = 0.95, na.rm = TRUE), 
      y > quantile(y, probs = 0.05, na.rm = TRUE), 
      y < quantile(y, probs = 0.95, na.rm = TRUE)
    )
  

  agg <-  ds$Canvas()$points(df, "x", "y")
  
  
  flip_180 <- function(x)  x[rev(seq_len(nrow(x))), ]
  invert <- function(x) 1 - x / max(x)
  
  agg$data %>% 
    flip_180() %>%
    log1p() %>%
    invert() %>%
    as.raster() %>% 
    plot()
  
})