Matt Owen

source code for "files/robots quadrature/nonuniform_quadrature.r"

return to portfolio
  1.  # matt owen
  2.  
  3.  
  4.  # chromosone object constructor
  5.  newchromosone <- function (v) {
  6.   obj <- c()
  7.   obj$genes <- sort(v)
  8.   attr(obj, 'class') <- 'chromosone'
  9.   obj
  10.  }
  11.  
  12.  # sampling functions # # # # # #
  13.  samp <- function (n, v, p) {
  14.   lis <- c()
  15.   v <- v/sum(v)
  16.   vec <- cumsum(v)
  17.   for (k in 1:n) {
  18.   vc <- vec
  19.   r <- runif(n=1, min=0, max=1)
  20.   vc <- (r < vc)*1:length(v)
  21.   vc <- vc[vc!=0]
  22.   lis <- c(lis, p[min(vc)])
  23.   }
  24.   lis
  25.  }
  26.  
  27.  sortselect <- function (p, f=function (x) { 1 }, n=1) {
  28.   fitness <- c()
  29.   for (n in 1:length(p)) {
  30.   fitness <- c(fitness, f(p[n]))
  31.   }
  32.   samp(n=n, fitness, p)
  33.  }
  34.  
  35.  # crossover functions # # # # # #
  36.  
  37.  # crossover average: average the two vectors
  38.  # should be a good mechanism for mixing solutions close to optimal
  39.  crossover_avg <- function (x, y) { newchromosone(sort((x$genes+y$genes)/2.)) }
  40.  
  41.  # crossover swap: swap two random entries
  42.  # adds variability to solutions
  43.  crossover_swap <- function (x, y) {
  44.   i <- floor(runif(n=1, min=0, max=length(x$genes)+1))
  45.   x$genes[i] <- y$genes[i]
  46.   newchromosone(sort(x$genes))
  47.  }
  48.  
  49.  
  50.  # mutate functions # # # # # #
  51.  
  52.  # mutate: naive implementation
  53.  mutate <- function (ch, mu=0.003) {
  54.   for (k in 1:length(ch$genes))
  55.   if (runif(n=1, min=0, max=1) < mu)
  56.   ch$genes[k] = runif(n=1, min=-5, max=5)
  57.   newchromosone(sort(ch$genes))
  58.  }
  59.  
  60.  # compute which partition a point belongs to
  61.  findpartition <- function (x, v) {
  62.   v <- sort(v)
  63.   maximum <- length(v)
  64.   v <- (x < v)*(1:length(v))-1
  65.   v <- v[v != -1]
  66.   min(v, maximum)+1
  67.  }
  68.  
  69.  # get a list of which segment something belongs to
  70.  findpartitions <- function (x, v) {
  71.   lis <- c()
  72.   for (el in x)
  73.   lis <- c(lis, findpartition(el, v))
  74.   lis
  75.  }
  76.  
  77.  
  78.  # compute which partition a point belongs to
  79.  findpartition <- function (x, v) {
  80.   v <- sort(v)
  81.   maximum <- length(v)
  82.   v <- (x < v)*(1:length(v))-1
  83.   v <- v[v != -1]
  84.   min(v, maximum)+1
  85.  }
  86.  
  87.  # get a list of which segment something belongs to
  88.  findpartitions <- function (x, v) {
  89.   lis <- c()
  90.   for (el in x)
  91.   lis <- c(lis, findpartition(el, v))
  92.   lis
  93.  }
  94.  
  95.  # constructor
  96.  newpiecewise <- function (f, p) {
  97.   p <- sort(p)
  98.   obj <- c()
  99.   obj$f <- c(f[[1]])
  100.   obj$p <- p
  101.   for (k in 1:length(p))
  102.   obj$f[[k+1]] <- f[[k+1]]
  103.   attr(obj, 'class') <- 'piecewisefunction'
  104.   obj
  105.  }
  106.  
  107.  # method to compute difference between functions
  108.  piecewisedifference <- function (obj, f, x) {
  109.   su <- 0
  110.   x <- sort(x)
  111.   y <- fapply(obj, x)
  112.   for (k in 1:(length(x)-1)) {
  113.   dx <- x[k+1] - x[k]
  114.   p2 <- abs(y[k+1] - f(x[k+1]))
  115.   p1 <- abs(y[k] - f(x[k]))
  116.   trap <- (p2+p1)/2*dx
  117.   su <- su+trap
  118.   }
  119.   su
  120.  }
  121.  
  122.  
  123.  # method to apply piecewise functions to points
  124.  fapply <- function (obj, x) {
  125.   partitions <- findpartitions(x, obj$p)
  126.   results <- c()
  127.   for (k in 1:length(x))
  128.   results <- c(results, obj$f[[partitions[[k]]]](x[k]))
  129.   results
  130.  }
  131.  
  132.  # linearize a function
  133.  linearize <- function (f, a, b) {
  134.   m <- (f(b) - f(a))/(b-a)
  135.   b <- f(a) - m * a
  136.   function (x) { m * x + b }
  137.  }
  138.  
  139.  # construct piecewise linear
  140.  piecewiselinear <- function (f, p) {
  141.   p <- sort(p)
  142.   flis <- c()
  143.   for (k in 1:(length(p)-1)) {
  144.   flis <- c(flis, linearize(f, p[k], p[k+1]))
  145.   }
  146.   return(newpiecewise(flis, p[2:(length(p)-1)]))
  147.  }
  148.  
  149.  # fitness
  150.  fit <- function (ch) {
  151.   y <- function (x) { x^3 }
  152.   fi <- piecewiselinear(y, c(-5, ch$genes, 5))
  153.   1/piecewisedifference(fi, y , seq(-5, 5, .1))
  154.  }
  155.  
  156.  
  157.  # program start # #
  158.  
  159.  # create the inital population
  160.  samples <- c()
  161.  x <- seq(-5, 5, .1)
  162.  y <- function (x) { x^3 }
  163.  
  164.  for (k in 1:10)
  165.   samples <- c(samples, newchromosone(sort(runif(n=4, min=-5, max=5))))
  166.  
  167.  lis <- c()
  168.  
  169.  val <- 0
  170.  m <- 0
  171.  best <- c()
  172.  
  173.  for (k in 1:length(samples)) {
  174.   lis <- c(lis, fit(samples[k]))
  175.   val <- fit(samples[k])
  176.   if (val > m) {
  177.   m <- val
  178.   best <- samples[k]
  179.   }
  180.  }
  181.  
  182.  pw <- piecewiselinear(function (x) { x^3 }, c(-5, best$genes, 5))
  183.  plot(x, fapply(pw, x), main="OLD BEST", xlab="x", ylab="y")
  184.  lines(x, y(x))
  185.  
  186.  
  187.  # start the loop
  188.  for (n in 1:20) {
  189.   # sample population
  190.   nextgeneration <- c()
  191.   parents <- sortselect(p=samples, f=fit, n=10)
  192.  
  193.   # mutate
  194.   #for (k in 1:length(parents))
  195.   # parents[k] <- mutate(parents[k])
  196.  
  197.   # mix
  198.   len <- length(parents)
  199.   for (k in 2:len) {
  200.   x <- parents[k-1]
  201.   y <- parents[k]
  202.   nextgeneration <- c(nextgeneration, crossover_avg(x, y))
  203.   }
  204.   nextgeneration <- c(nextgeneration, crossover_avg(parents[1], parents[length(parents)]))
  205.   samples <- nextgeneration
  206.  }
  207.  
  208.  for (k in 1:length(samples)) {
  209.   lis <- c(lis, fit(samples[k]))
  210.   val <- fit(samples[k])
  211.   if (val > m) {
  212.   m <- val
  213.   best <- samples[k]
  214.   }
  215.  }
  216.  
  217.  # text output
  218.  m <- 0
  219.  best <- c()
  220.  val <- 0
  221.  for (k in 1:length(samples)) {
  222.   lis <- c(lis, fit(samples[k]))
  223.   val <- fit(samples[k])
  224.   if (val > m) {
  225.   m <- val
  226.   best <- samples[k]
  227.   }
  228.  }
  229.  
  230.  pw <- piecewiselinear(function (x) { x^3 }, c(-5, best$genes, 5))
  231.  x <- seq(-5, 5, .1)
  232.  y <- function (x) { x^3 }
  233.  plot(x, fapply(pw, x), main="NEW BEST", xlab="x", ylab="y")
  234.  lines(x, y(x))
  235.  
  236.  
  237.  
  238.  # output image file
  239.  
  240.  
  241.  #findpartitions(c(-1.5, -.5, .5, 1.5), c(-1, 0, 1))