R でエンタメ

---

数独†

• Rで数独

カードゲーム†

• gamesNws: Playing games using a NWS Server
• A pure R poker hand evaluator
• Playing cards, with R
• Playing cards in Vegas?
• Simulation of Blackjack: the odds are not with you

迷路作成†

迷路を抜けろ！

maze <- function(xsize=150, ysize=150, color="black") # 河童の屁
{
  old <- par(mar=c(1,1,1,1))
  xsize3 <- xsize+3
  ysize3 <- ysize+3
  plot(c(2, xsize3), c(2, ysize3), type="n", axes=FALSE, xlab="", ylab="", asp=1)
  p <- matrix(FALSE, ysize+4, xsize+4)
  p[3:(ysize+2), 3:(xsize+2)] <- TRUE
  segments(c(2, 2, 3, xsize3), c(2, ysize3, 2, 2),
         c(2, xsize3-1, xsize3, xsize3), c(ysize3, ysize3, 2, ysize3), col=color)
  n <- xsize*ysize
  n0 <- n%/%20
  istart <- 2
  check.all <- TRUE
  repeat {
    if(n > n0) {
      repeat {
        i1 <- sample(2:ysize3, 1)
        j1 <- sample(2:xsize3, 1)
        if (p[i1, j1] == FALSE && (i1 != 2 || j1 != 2) && (i1 != ysize3 || j1 != xsize3)) {
          break
        }
      }
    }
    else {
      found <- FALSE
        for (i1 in istart:ysize3) {
        for (j1 in 2:xsize3) {
          if (p[i1+1, j1] || p[i1, j1+1]) {
            found <- TRUE
            break
          }
        }
        if(found) {
          istart <- i1
          break
        }
      }
      check.all <- FALSE
    }
    repeat {
      ix <- iy <- integer(3)
      m <- 0
      if (p[i1+1, j1]) {
        m <- m+1
        ix[m] <- 1
      }
      if (p[i1, j1+1]) {
        m <- m+1
        iy[m] <- 1
      }
      if (check.all) {
        if (p[i1-1, j1]) {
          m <- m+1
          ix[m] <- -1
        }
        if (p[i1, j1-1]) {
          m <- m+1
          iy[m] <- -1
        }
      }
      if (m == 0) {
        break
      }
      else if (m >= 2) {
        m <- sample(m, 1)
      }
      i2 <- i1+ix[m]
      j2 <- j1+iy[m]
      segments(j1, i1, j2, i2, col=color)
      p[i2, j2] <- FALSE      
      i1 <- i2
      j1 <- j2
      n <- n-1
    }
    if(n == 0) {
      break
    }
  }
  par(old)
}

N Queen†

N x N チェスボードに，N 個のクイーンを平和共存させよ！

N.Queen <- function(N=8) # 河童の屁
{
  Solve <- function(i)
  {
    if (i > N) {
      Counter <<- Counter+1
      cat(sprintf("\nNo. %i\n", Counter))
      Board <- matrix(".", N, N)
      Board[cbind(1:N, Row)] <- "Q"
      colnames(Board) <- letters[1:N]
      rownames(Board) <- N:1
      print(Board, quote=FALSE)
      return
    }
    for (j in 1:N) {
      if (Col[j] == TRUE && RightUp.LeftDown[i+j-1] == TRUE && RightDown.LeftUp[i-j+N] == TRUE ) {
        Row[i] <<- j
        Col[j] <<- RightUp.LeftDown[i+j-1] <<- RightDown.LeftUp[i-j+N] <<- FALSE
        Solve(i+1)
        Row[i] <<- Col[j] <<- RightUp.LeftDown[i+j-1] <<- RightDown.LeftUp[i-j+N] <<- TRUE
      }
    }
  }
  Row <- integer(N)
  Col <- rep(TRUE, N)
  RightUp.LeftDown <- RightDown.LeftUp <- rep(TRUE, 2*N-1)
  Counter <- 0
  Solve(1)
}

N.Queen(12) の解の最後のもの

No. 14200
   a b c d e f g h i j k l
12 . . . . . . . . . . . Q
11 . . . . . . . . . Q . .
10 . . . . . . . Q . . . .
9  . . . . Q . . . . . . .
8  . . Q . . . . . . . . .
7  Q . . . . . . . . . . .
6  . . . . . . Q . . . . .
5  . Q . . . . . . . . . .
4  . . . . . . . . . . Q .
3  . . . . . Q . . . . . .
2  . . . Q . . . . . . . .
1  . . . . . . . . Q . . .

ハノイの塔†

n 枚の板を a から b へ，作業場所として c を使って移せ！

Hanoi <- function(n=15) # 河童の屁
{
  move <- function(a, b, c, n)
  {
    if (n) {
      move(a, c, b, n-1)
      animation(a, b)
      move(c, b, a, n-1)
    }
  }

  animation <- function(a, b)
  {
  #  cat("from", a, "to", b, "\n")
    for (i in 1:n) {
      if (hanoi[i, a] != 0) {
        disc <- hanoi[i, a]
        hanoi[i, a] <<- 0
        draw.disc(a, i, disc, FALSE)
        break
      }
    }
    for (i in n:1) {
      if (hanoi[i, b] == 0) {
        hanoi[i, b] <<- disc
        draw.disc(b, i, disc, TRUE)
        break
      }
    }
  }

  draw.disc <- function(a, i, disc, flag)
  {
    x <- n+2*n*(a-1)
    y <- n+1-i
    rect(x-disc, y-0.45, x+disc, y+0.45, col=if (flag) "blue" else "white", border="white")
  }

  hanoi <- cbind(1:n, rep(0, n), rep(0, n))
  old <- par(mar=c(1,1,1,1))
  plot(c(0, 6*n), c(0, n+1), type="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
  text(c(n, 3*n, 5*n), -0.5, letters[1:3], xpd=TRUE)
  for (i in 1:n) {
    draw.disc(1, i, i, TRUE)
  }
  move(1, 2, 3, n)
  par(old)
}
Hanoi()

作業の途中

R の唄†

• The Most Romantic Electro-Grunge Statistical Computing Song Ever Made

RLastFM (R interface to last.fm API)†

• RLastFM (R interface to last.fm API)

Rでスポーツ統計†

バットマン方程式- バットサインの数式による描画†

the batman equation ggplot2 利用†

plot 関数で†

ggplot2 版の 40 倍の速度 by 河童の屁は，河童にあらず，屁である。

plot(0, 0, xlim=c(-7, 7), ylim=c(-3, 3), type="n", bty="n",
     axes=FALSE, xlab="", ylab="", asp=1)
x11 <- seq(3, 7, 0.001)
x12 <- seq(7, 4, -0.001)
y1 <- c(3*sqrt(1-(x11/7)^2), -3*sqrt(1-(x12/7)^2))
x2 <- seq(4, 0, -0.001)
y2 <- x2/2-(3*sqrt(33)-7)*x2^2/112-3 + sqrt(1-(abs(x2-2)-1)^2)
x345 <- c(1, 0.75, 0.5, -0.5, -0.75, -1, 1)
y345 <- c(1, 3, 2.25, 2.25, 3, 1, 1)
x6 <- seq(1, 3, 0.001)
y6 <- 1.5-0.5*x6+6*sqrt(10)*(1/7-sqrt(4-(x6-1)^2)/14)
lines(c(x345, x6, x11, x12, x2, NA, -x2, NA, -x11, -x12, NA, -x6),
      c(y345, y6, y1,       y2, NA,  y2, NA,  y1,        NA,  y6))

やはり，塗りつぶされていないとね。

Le Mond の数学パズルをRで解く†

• Le Monde puzzle [1]

エレガントな解を探しましょう†

規則に従う要素を持つ正方行列（出題：河童の屁は，河童にあらず，屁である。）†

n = 2 のとき
1 2
4 3
n = 3 のとき
1 2 3
8 9 4
7 6 5
n = 4 のとき
 1  2  3  4
12 13 14  5
11 16 15  6
10  9  8  7

のように，時計回りの渦のように 1:(n^2) を配置する配列を返す関数 a <- function(n) { } の中身を書く。
どんなライブラリ，関数を使っても良い。
空白類（空白，タブ，改行）を除く文字数が最小のものを最良のプログラムとする。

まあ，面白くもないということのようですが，例解?。174 文字。2011/09/10

a <- function(n)
{
	x <- diag(n)
	for (i in 1:((n+1)%/%2)) {
		s <- n+1-2*i
		if (s == 0) {
			x[i, i] <- n^2
			break
		}
		j <- rep(i, s)
		k <- c(i:(i+s-1), j+s, i+s:1)
		x[cbind(c(j, k), c(k, j))] <- n^2-(s+1)^2+1:(4*s)
	}
	return(x)
}

これはどうでしょう？　143 文字。2011/10/26

b <- function(n)
{
        A <- function(m){
                x <- ncol(m)
		  y <- 2*x+1
                rbind(1:(x+1), cbind(m[x:1,x:1] + y, (x+2):y))
        }
        B <- function(m){ 
                if(ncol(m) < n) B(A(m))
                else m
        }
        B(diag(1))
}

これは156文字だけど、もう少し短くできるかな？

b <- function(n)
{
        A <- function(m){
                x <- ncol(m)
		  y <- 2*x+1
                rbind(1:(x+1), cbind(m[x:1,x:1] + y, (x+2):y))
        }
	eval( parse(text = paste(rep(c('A(','diag(1)',')'),c(n-1,1,n-1) )))) 
}

あるいはシンプルに for loop で・・・　136文字。

b <- function(n)
{
        A <- function(m){
                x <- ncol(m)
		  y <- 2*x+1
                rbind(1:(x+1), cbind(m[x:1,x:1] + y, (x+2):y))
        }
	out <- diag(1)
	if(n>1) for(i in 2:n) out <- A(out)
	out
}

文字数最小なのだから，変数名も1文字にすれば，この解は 128 文字解ですね。

ははは。うっかりしていました。アウトですね。　<- を = にするというのはなんとなく自分で許せませんでした。

お題：文字列を先頭から見て同じところまで除去†

func <- function(s) # 河童の屁
{
	x <- lapply(s, function(x) unlist(strsplit(x, "")))
	y <- sapply(x, "[",  1:min(sapply(x, length)))
	z <- sum(cumsum(apply(y, 1, function(y) length(unique(y)) != 1)) == 0)
	sapply(s, function(x) sub(sprintf(".{%i}", z), "", x))
}

実行結果

> func(c("abc123", "abcdefg", "abrt"))
 abc123 abcdefg    abrt 
 "c123" "cdefg"    "rt" 
> func(c("abc123", "abcdefg", "wabrt"))
   abc123   abcdefg     wabrt 
 "abc123" "abcdefg"   "wabrt" 
> func(c("かみさま", "かみゆい", "かみさん", "かみかぜ"))
かみさま かみゆい かみさん かみかぜ 
  "さま"   "ゆい"   "さん"   "かぜ" 
> func(c("かみさま", "かみゆい", "かみさん", "abc"))
  かみさま   かみゆい   かみさん        abc 
"かみさま" "かみゆい" "かみさん"      "abc"

お題：アラビア数字・ローマ数字変換†

R には as.roman があるけど，変換例に挙げられているのをやってみると as.roman(3999) は NA になる。これじゃまずいなあ。また，逆変換はない。
誰かやってみる？）
プロジェクトオイラーの問題89も参照。

• まあ，余りエレガントでもないけど。

  encode <- function(r) # ローマ数字へ
  {
  	str <- NULL
  	s <- as.integer(rev(c(rep("0", 3), unlist(strsplit(as.character(r), ""))))[1:4])
  	s4 <- s[4]
  	s3 <- s[3]
  	s2 <- s[2]
  	s1 <- s[1]
  	if (s4 > 0) str <- c(str, rep("M", s4))
  
  	if (s3 == 9) str <- c(str, "C", "M")
  	else if (s3 > 5) str <- c(str, "D", rep("C", s3-5))
  	else if (s3 == 5) str <- c(str, "D")
  	else if (s3 == 4) str <- c(str, "C", "D")
  	else if (s3 > 0) str <- c(str, rep("C", s3))
  
  	if (s2 == 9) str <- c(str, "X", "C")
  	else if (s2 > 5) str <- c(str, "L", rep("X", s2-5))
  	else if (s2 == 5) str <- c(str, "L")
  	else if (s2 == 4) str <- c(str, "X", "L")
  	else if (s2 > 0) str <- c(str, rep("X", s2))
  
  	if (s1 == 9) str <- c(str, "I", "X")
  	else if (s1 > 5) str <- c(str, "V", rep("I", s1-5))
  	else if (s1 == 5) str <- c(str, "V")
  	else if (s1 == 4) str <- c(str, "I", "V")
  	else if (s1 > 0) str <- c(str, rep("I", s1))
  	return(paste(str, collapse=""))
  }
  decode <- function(r) # ローマ数字から
  {
  	v <- 0
  	str <- unlist(strsplit(r, ""))
  	repeat {
  		n <- length(str)
  		if (n == 0) break
  		s <- str[1]
  		if (grepl(s, "IXC") && n >= 2) {
  			s2 <- str[2]
  			if (s == "I" && s2 == "V") {
  				v <- v+4
  				str <- str[-(1:2)]
  				next
  			}
  			else if (s == "I" && s2 == "X") {
  				v <- v+9
  				str <- str[-(1:2)]
  				next
  			}
  			else if (s == "X" && s2 == "L") {
  				v <- v+40
  				str <- str[-(1:2)]
  				next
  			}
  			else if (s == "X" && s2 == "C") {
  				v <- v+90
  				str <- str[-(1:2)]
  				next
  			}
  			else if (s == "C" && s2 == "D") {
  				v <- v+400
  				str <- str[-(1:2)]
  				next
  			}
  			else if (s == "C" && s2 == "M") {
  				v <- v+900
  				str <- str[-(1:2)]
  				next
  			}
  		}
  		if (s == "I") {
  			v <- v+1
  		}
  		else if (s == "V") {
  			v <- v+5
  		}
  		else if (s == "X") {
  			v <- v+10
  		}
  		else if (s == "L") {
  			v <- v+50
  		}
  		else if (s == "C") {
  			v <- v+100
  		}
  		else if (s == "D") {
  			v <- v+500
  		}
  		else if (s == "M") {
  			v <- v+1000
  		}
  		str <- str[-1]
  	}
  	return(v)
  }
  実行例
  > as.roman(3999)
  [1] <NA>
  > encode(3999)
  [1] "MMMCMXCIX"
  > decode("MMMCMXCIX")
  [1] 3999

錯視のカタログ†

北岡明佳の錯視のページ 色々な例を見ているだけで面白い
Rで作図してみる例題としては，高度のテクニックが必要のものもあるようだ（何で作図しているのだろうかなあ，Rではなさそう(笑)

• ↑　これはRと何の関係が...。誰かコード書いてってこと？
• 関係ないと思う人はそう思えばよい。
  プログラムを書くのがエンターテインメントの人もいるんじゃないかなと思うけど。
  Rで描いたらどうなるかという，練習問題としてとらえれば良いでしょう。ここは「エンタメ」なんだから，いろいろなエンタメがあって良いと思うけど？
  「R を使ってなんとかかんとか」というのがありなら，「何とかかんとかを R で描く（書く）」というのがあって何の不都合があるんでしょうか？頼まれてもかけない人は黙っていればよい。
• ちょっと書いてみました。

  factor <- 0.04
  n <- 20
  t1 <- seq(0, 2*pi, length=13)
  t3 <- seq(0, 2*pi, length=801)[1:n]
  #plot(l2*cos(t2), l2*sin(t2))
  color <- c("black", "white")
  color2 <- c("#cfc947", "#4e79f5")
  color2 <- c("aquamarine4", "green1")
  color2 <- c("aquamarine4", "gold2")
  
  parts <- function(x0, y0, i)
  {
  	l <- 0.8^(i-1)
  	j <- i%%2+1
  	l4 <- rep(l, n)
  	l4 <- c(l4, l4-l*0.8)
  	for (t in seq(0, 2*pi, length=41)[-1]) {
  		t4 <- t+c(t3, rev(t3))
  		polygon(l4*cos(t4)+x0, l4*sin(t4)+y0, border=color[j], col=color[j])
  		j <- 3-j
  	}
  	t2 <- atan(l*factor*sin(t1)/(0.8*l+0.1*l*cos(t1)))
  	l2 <- sqrt((0.9*l+0.1*l*cos(t1))^2+(factor*sin(t1))^2)
  	for (t in seq(0, 2*pi, length=41)[-1]) {
  		polygon(l2*cos(t+t2)+x0, l2*sin(t+t2)+y0, border=color2[j], col=color2[j])
  		j <- 3-j
  	}
  }
  old <- par(mar=rep(0, 4))
  plot(0, 0, type="n", xlim=c(0, 8), ylim=c(0, 6), asp=1, axes=FALSE, xlab="", ylab="")
  draw <- function(x0, y0)
  {
  	for (i in 1:15) {
  		parts(x0, y0, i)
  	}
  }
  for (x0 in 1:4*2-1) for (y0 in 1:3*2-1) draw(x0, y0)
  for (x0 in 1:3*2) for (y0 in 1:2*2) draw(x0, y0)
  text(5, -0.4, "「蛇の回転」Copyright Akiyoshi Kitaoka 2003 (September 2, 2003)", pos=3, xpd=TRUE)
  par(old)

ハロウィン†

• R でハロウィンカード [#ye3ef3e8]

チェス†

• First attempt at Chess Data Mining

ロゴ作成†

• Using R to create a logo: Simple

くじ引き†

• New Powerball (lottery) Rules Will Cost You More
• 宝くじがどのくらい当たるか実感してみよう

Pasta Geometrics†

• Pasta Geometrics こんなのでよいのかなあ? 河童の屁は，河童にあらず，屁である。
  Agnolotti

  PI <- THETA <- KAPPA <- matrix(0, 61, 101)
  for (i in 0:60) {
      for (j in 0:100) {
          PI[i, j] <- (10 * sin(i * pi / 120) ^ 0.5 + i / 400 * sin(3 * j * pi / 10)) *
              cos(19 * j * pi / 2000 + 0.03 * pi)
          THETA[i, j] <- (10 * sin(i * pi / 120) + i / 400 * cos(3 * j * pi / 10)) *
              sin(19 * j * pi / 2000 + 0.03 * pi)
          KAPPA[i, j] <- 5 * cos(i * pi / 120) ^ 5 * sin(j * pi / 100) -
              5 * sin(j * pi / 100) * cos(i * pi / 120) ^ 200
      }
  }
  library(rgl)
  plot3d(c(PI), c(THETA), c(KAPPA), aspect = c(1, 1, 0.25)) 

  

plot3d は隠れた部分の消去機能がないので，見た目が悪い。自分で回転させれば形がよくわかる。
Fagottini

PI <- THETA <- KAPPA <- matrix(0, 201, 51)
for (i in 0:200) {
    for (j in 0:50) {
        alpha <- (0.8 + sin(i * pi / 100) ^ 8 - 0.8 * cos(i * pi / 25)) ^ 1.5 + 0.2 + 0.2 * 
            sin(i * pi / 100)
        beta <- (0.9 + cos(i * pi / 100) ^ 8 - 0.9 * cos((4 * i + 3) * pi / 100)) ^ 1.5 + 
            0.3 * cos(i * pi / 100)
        gamma <- 4 - 4 * j / 500 * (1 + cos(i * pi / 100) ^ 8 - 0.8 * cos(i * pi / 25)) ^ 1.5
        PI[i, j] <- cos(i * pi / 100) * (alpha * sin(j * pi / 100) ^ 8 +  
            0.6 * (2 + sin(i * pi / 100) ^ 2) * sin(j * pi / 50) ^ 2)
        THETA[i, j] <- sin(i * pi / 100) * (beta * sin(j * pi / 100) ^ 8 +  
            0.6 * (2 + cos(i * pi / 100) ^ 2) * sin(j * pi / 50) ^ 2)
        KAPPA[i, j] <- (1 + sin((j - 50) * pi / 100)) * 
            (gamma - 4 * j / 500 * (1 + sin(i * pi / 100) ^ 8 - 0.8 * cos(i * pi / 25)) ^ 1.5)
    }
}
library(rgl)
plot3d(c(PI), c(THETA), c(KAPPA), aspect = c(1, 1, 0.5)) 

CX's 'reality distortion' pie-graph (フジTV のステマ円グラフ)†

http://blog.livedoor.jp/dqnplus/archives/1692812.html

これは、ひどい。 が、遊びで描いてみた。

library(grid)
png(filename = "~/Desktop/fuji-pie.png", width = 360, height = 360, units = "px")
num <- c(94, 97, 78, 78)
distortion <- c(1, 1.2, 1, 1)
dist_num <- num * distortion
cum <- c(0, cumsum(dist_num))
rad <- cum/sum(dist_num) * 360/180 * pi
colors <- c("#0000FFAA", "#FF0000AA", "#000055AA", "#0000AAAA")
labels <- paste(c("50ies\n", "10~20ies\n", "30ies\n", "40ies\n"), 
    num, " ps", sep = "")
lab_cex <- c(1, 1.5, 1, 1)
lab_col <- c("white", "yellow", "white", "white")
div <- 200  # 円弧部分の分割数
dx <- 0     # 半径=0.5 として、x 方向への中心のズレ
dy <- 0.1   # 半径=0.5 として、y 方向への中心のズレ
for (i in 1:length(num)) {
    x <- c(0.5 + dx, 0.5 + 0.5 * sin(rad[i]), 0.5 + 0.5 * sin(seq(rad[i], 
        rad[i + 1], length.out = div)), 0.5 + 0.5 * sin(rad[i + 1]))
    y <- c(0.5 + dy, 0.5 + 0.5 * cos(rad[i]), 0.5 + 0.5 * cos(seq(rad[i], 
        rad[i + 1], length.out = div)), 0.5 + 0.5 * cos(rad[i + 1]))
    grid.polygon(x = x, y = y, gp = gpar(fill = colors[i]))
    grid.text(x = 0.5 + 0.4 * sin((rad[i] + rad[i + 1])/2), y = 0.5 + 0.4 * 
        cos((rad[i] + rad[i + 1])/2), label = labels[i], gp = gpar(col = lab_col[i], 
        cex = lab_cex[i]))
}
dev.off() 

古典的グラフィックでも。

fuji.pie <- function(x, center, labels, colors, text.colors, text.cex) {
    if (missing(center)) 
        center <- c(0, 0.3)
    i <- (1:360)/360 * 2 * pi
    arc <- cut(1:360, c(0, round(cumsum(x) * 360/sum(x)), Inf))
    plot.new()
    plot.window(xlim = c(-1, 1), ylim = c(-1, 1), "", asp = 1)
    l <- levels(arc)
    for (j in 1:length(x)) {
        k <- c(i[arc == l[j]], max(i[arc == l[j]]) + 1/360 * 2 * pi)
        polygon(c(center[1], cos(k)), c(center[2], sin(k)), col = colors[j])
        text(mean(cos(k))/1.5, mean(sin(k))/1.5, labels[j], cex = text.cex[j], 
            col = text.colors[j])
    }
}
fuji.pie(c(94, 78, 78, 97), labels = c("50代\n94人", "40代\n78人", 
    "30代\n78人", "10〜20代\n97人"), colors = c("deepskyblue", "dodgerblue1", 
    "dodgerblue3", "brown1"), text.colors = c("white", "white", "white", 
    "yellow"), text.cex = c(1.5, 1.5, 1.5, 3)) 

• >>fuji.pie() お見事、完敗です。でも、この関数を誰かが使うのか?? (笑)　エンタメだからいいですよね。

Music†

• Music player in R (Linux) [#v642664e]
• Programming instrumental music from scratch [#wd100e3a]

ウォーリーを探せ†

• [raster パッケージによるもの:http://stackoverflow.com/questions/8563604/how-to-find-waldo-with-r]]
• Rでウォーリーを探してみた
• Where's Waldo? Image Analysis in R

Hangman†

• Hangman in R: A learning experience

Minesweeper†

• Minesweeper

カッシーニの卵形線†

• Animation basics for a vacation

スピログラフ†

• Spirograph with R
• スピログラフ

Solving 9-puzzle with GNU R†

星空†

pchでいつの間にか全角文字が使えるようになっていたのを知ったので、お昼休みに「★」で星空を作成してみた -- 谷村 2013-02-22 (金) 12:48:52

• 星の位置は、CSR (完全空間ランダム)による分布
• 星の大きさは対数正規分布
• 星の色はここを参照。でも、結果的に適当 :-p
• filled.contour()は、keyを内部でハードコーディングしているので、オプションでkeyを取り除けず、仕方がないのでpngパッケージで切り落とした

library(splancs)
library(png)
f <- paste0(tempfile(), ".png")
png(file = f)
op <- par(bg = "#000022", mai = c(0, 0, 0, 0))
filled.contour(volcano, nlevels = 100,
               color.palette = colorRampPalette(c("#000000", "#007BC3"),
                 space = "Lab", bias = .5))
n <- 1000
p <- matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0), ncol = 2, byrow = TRUE)
cols <- c(rep("#FFFFFF", 10), "#FAFAFA", "#FFEF85", "#FFB74C")
points(csr(p, n), pch = "★", col = cols, cex = rlnorm(n, 0, .5) * .2)
par(op)
dev.off()
p <- readPNG(f)
unlink(f)
p1 <- p[, -(440:480), ]
writePNG(p1, "20130222.png")

Happy St Patrick’s Day†

Rでスネークゲーム†

Rで切符問題†

 rv.kippu<-function(num){
 exp <- NULL
 stackpointer <- 0
 stack <- list()
 
 perm<-function(v){
   temp<-NULL
   if(length(v) == 1){
     v[1]
   }else{
     for(i in 1:length(v)){
       temp <- rbind(temp,cbind(v[i],Recall(v[-i])))
       
     }
     temp
   }
   
 }  
 push <- function(d){
   stack[[stackpointer + 1]] <<- d
   stackpointer <<- (stackpointer + 1)
   
 }
 pop <- function(){
   if(stackpointer == 0){
     NULL
   }else{
     value <- stack[[stackpointer]]
     stack[[stackpointer]] <<- NULL
     stackpointer <<- (stackpointer - 1)
     value
     }
 }  
 rpn2normal<-function(v){
   len<-length(v)
   for(i in 1:len){
     if(any(grep("^[-+*/]",v[i])) ){
       n2<-pop()
       n1<-pop()
       exp<-paste("(",n1,v[i],n2,")",sep="")
       push(exp)
     }else{
       if(any(grep("[0-9]",v[i]))) {
         push(v[i])
       }
     }
   }
   invisible(exp)
 }
 if(num >= 0 && num <= 9999){
   n1 <- as.character(n1 <- num %/% 1000)
   n2 <- as.character(n2 <- num %% 1000 %/% 100)
   n3 <- as.character(n3 <- num %% 100 %/% 10)
   n4 <- as.character(n4 <- num %% 10)
   t<-perm(c(n1,n2,n3,n4))
   op <- c("+","-","*","/")
   exps<-NULL  
   # num num opr num num opr opr
   for(i in c(1:24)){
     nums<-t[i,]
     tmp<-NULL
     for(i1 in op){
       for(i2 in op){
         for(i3 in op){
           tmp<-c(nums[1],nums[2],i1,nums[3],nums[4],i2,i3)
           exps<-c(exps,tmp)
         }
       }
     }
   }
   # num num opr num opr num opr
   for(i in c(1:24)){
     nums<-t[i,]
     tmp<-NULL
     for(i1 in op){
       for(i2 in op){
         for(i3 in op){
           tmp<-c(nums[1],nums[2],i1,nums[3],i2,nums[4],i3)
           exps<-c(exps,tmp)
         }
       }
     }
   }
   # num num num opr  num opr opr
   for(i in c(1:24)){
     nums<-t[i,]
     tmp<-NULL
     for(i1 in op){
       for(i2 in op){
         for(i3 in op){
           tmp<-c(nums[1],nums[2],nums[3],i1,nums[4],i2,i3) 
           exps<-c(exps,tmp)
         }
       }
     }
   }
   # num num num opr num opr opr     
   for(i in c(1:24)){
     nums<-t[i,]
     tmp<-NULL
     for(i1 in op){
       for(i2 in op){
         for(i3 in op){
           tmp<-c(nums[1],nums[2],nums[3],i1,nums[4],i2,i3) 
           exps<-c(exps,tmp)
         }
       }
     }
   }
   # num num opr num num opr opr 
   for(i in c(1:24)){
     nums<-t[i,]
     tmp<-NULL
     for(i1 in op){
       for(i2 in op){
         for(i3 in op){
           tmp<-c(nums[1],nums[2],i1,nums[3],nums[4],i2,i3) 
           exps<-c(exps,tmp)
         }
       }
     }
   }
  tmp <- apply(matrix(exps,byrow=T,ncol=7),1,rpn2normal)
   ret <- NULL
   for(i in tmp){
     #print(class(i))
     ans <- eval(parse(text=i))
     
     if(is.finite(ans) &&(ans == 10)){
       s <- paste(i,"=",ans)
       #print(s)
       ret <- c(ret,s)
     }
  }
   ret <- matrix(ret,ncol=1)
   return(unique(ret))
 }else{
   print("please input 4 digits' number" )    
 }}# by Mamekkoro

• if(any(grep("^[-\+\*/]",v[i])) ){ は if(any(grep("^[-+*/]",v[i])) ){ でなきゃ，エラーになる。
  1111 など，不可能な場合に，「以下にエラー matrix(ret, ncol = 1) : 'data' はベクトル型でなくてはなりませんが、'NULL' でした 」エラーを吐くのはみっともない。
• 了解しました by mamekkoro
• 定型的な処理を繰り返すのに for を多用するのも善し悪し。少なくとも，i1, i2, i3 の op は，3 重の for ループに共通する 4^3 × 3 の matrix で表すことができる。
  計算時間はたいして違わないが，視認性には大きな違いがある。
  2.5倍くらい速くなった。
  簡潔にしたら，いろいろ不要な部分や間違い？も発見できた。目を通す部分は少ないに越したことはない。

  rv.kippu2 <- function(num) {
      perm <- function(v) {
          temp <- NULL
          if (length(v) == 1) {
              v[1]
          } else {
              for (i in 1:length(v)) {
                  temp <- rbind(temp, cbind(v[i], Recall(v[-i])))
              }
              temp
          }
      }
      rpn2normal <- function(v) {
          stackpointer <- 0
          stack <- character(10)
          for (i in seq_along(v)) {
              if (grepl("[-+*/]", v[i])) {
                  n2 <- stack[stackpointer]
                  stackpointer <- stackpointer - 1
                  n1 <- stack[stackpointer]
                  exp <- paste("(", n1, v[i], n2, ")", sep = "")
                  stack[stackpointer] <- exp
              } else if (grepl("[0-9]", v[i])) {
                  stackpointer <- stackpointer + 1
                  stack[stackpointer] <- v[i]
              }
          }
          invisible(exp)
      }
      if (num < 1000 || num > 9999) 
          stop("please input 4 digits' number")
      temp <- unlist(strsplit(as.character(num), ""))
      temp <- perm(temp)
      nummat <- apply(temp, 2, rep, each = 64)
      op <- c("+", "-", "*", "/")
      opmat <- sapply(expand.grid(op, op, op), as.character)
      opmat <- apply(opmat, 2, rep, 24)
      exps5 <- cbind(nummat, opmat)            # num num num num opr opr opr
      exps1 <- exps5[, c(1:2, 5, 3:4, 6:7)]    # num num opr num num opr opr
      exps2 <- exps5[, c(1:2, 5, 3, 6, 4, 7)]  # num num opr num opr num opr
      exps3 <- exps5[, c(1:3, 5:6, 4, 7)]      # num num num opr opr num opr
      exps4 <- exps5[, c(1:3, 5, 4, 6:7)]      # num num num opr num opr opr
      exps <- rbind(exps1, exps2, exps3, exps4, exps5)
      ans <- apply(exps, 1, rpn2normal)
      ans <- sapply(ans, function(e) eval(parse(text = e)))
      ans <- ans[ans == 10]
      return(unique(names(ans)))
  }  # by Mamekkoro. modified by  Kappa-no-He

  さらに 1.8 倍ほど速くなった。

  rv.kippu3 <- function(num) {
      perm <- function(v) {
          temp <- NULL
          if (length(v) == 1) {
              v[1]
          } else {
              for (i in 1:length(v)) {
                  temp <- rbind(temp, cbind(v[i], Recall(v[-i])))
              }
              temp
          }
      }
      rpn2normal <- function(v) {
          return(rbind(
          sprintf("%s%s(%s%s(%s%s%s))", v[1], v[5], v[2], v[6], v[3], v[7], v[4]),
          sprintf("(%s%s%s)%s(%s%s%s)", v[1], v[5], v[2], v[6], v[3], v[7], v[4]),
          sprintf("((%s%s%s)%s%s)%s%s", v[1], v[5], v[2], v[6], v[3], v[7], v[4]),
          sprintf("(%s%s(%s%s%s))%s%s", v[1], v[5], v[2], v[6], v[3], v[7], v[4]),
          sprintf("%s%s((%s%s%s)%s%s)", v[1], v[5], v[2], v[6], v[3], v[7], v[4])))
      }
      if (num < 1000 || num > 9999) 
          stop("please input 4 digits' number")
      temp <- unlist(strsplit(as.character(num), ""))
      temp <- perm(temp)
      nummat <- apply(temp, 2, rep, each = 64)
      op <- c("+", "-", "*", "/")
      opmat <- sapply(expand.grid(op, op, op), as.character)
      opmat <- apply(opmat, 2, rep, 24)
      exps <- cbind(nummat, opmat)
      ans <- apply(exps, 1, rpn2normal)
      ans <- sapply(ans, function(e) eval(parse(text = e)))
      ans <- ans[!is.na(ans) & ans == 10]
      return(unique(names(ans)))
  }  # by Mamekkoro. modified by  Kappa-no-He

• すごい簡潔になりました。勉強になります。ありがとうございました。次はHitandBlowを作ってみようともいます。mamekkoro

映画†

• Predicting movie ratings with IMDb data and R
• Who knows the Oscar winners? The betting markets, probably.
• Top 250 Movies at IMDb

覆面算†

smm <- function(str) {
    str <- tolower(gsub("[ \t]", "", str))
    str <- unlist(strsplit(str, ""))
    if (grepl(str[1], "+-")) {
        stop("Unary minus or plus is invalid.")
    }
    str <- c(str, ".")
    nChar <- length(str)
    operatorPosition <- NULL
    for (i in seq_len(nChar)) {
        if (match(str[i], c("+", "-", "=", "."), 0) != 0) {
            operatorPosition <- c(operatorPosition, i)
        } else if (match(str[i], letters, 0) == 0) {
            stop("Invalid character(s).")
        }
    }
    nTerm <- length(operatorPosition)
    word <- vector("list", nTerm)
    op <- rep("+", nTerm + 1)
    start <- 1
    for (i in seq_len(nTerm)) {
        position <- operatorPosition[i]
        word[[i]] <- str[start:(position - 1)]
        op[i + 1] <- str[position]
        start <- position + 1
    }
    charSet <- unique(unlist(word))
    if (length(charSet) > 10) {
        stop("Too many characters.")
    }
    library(e1071)
    perm <- t(e1071::permutations(10) - 1)
    term <- matrix(0, nTerm, ncol(perm))
    nonZeroTop <- rep(TRUE, ncol(perm))
    for (i in seq_len(nTerm)) {
        subscript <- match(word[[i]], charSet)
        term[i, ] <- colSums(10^(length(subscript):1-1) * perm[subscript, ])
        if (op[i] == "-") {
            term[i, ] <- -term[i, ]
        }
        nonZeroTop <- nonZeroTop & perm[subscript[1], ] != 0
    }
    leftsideTerm <- colSums(term[1:(nTerm - 1), ])
    answerFlag <- leftsideTerm == term[nTerm, ] & nonZeroTop
    ans <- unique(t(term[, answerFlag]))
    if (nrow(ans) == 0) {
        "No solution."
    } else if (nrow(ans) == 1) {
        op <- op[-nTerm-1]
        op[1] <- ""
        cat(gsub("--", "-", paste(op, ans, collapse="", sep="")))
    } else {
        ans
    }
}

実行例

> smm("Send + More = Money")
9567+1085=10652
> smm("Money - Send = More")
10652-9567=1085
> smm("Donald+Gerald=Robert")
526485+197485=723970
> smm("Bill+William+Monica=Clinton")
9600+1600634+457623=2067857
> smm("Green+Orange=Colors")
83446+135684=219130
> smm("Manet+Matisse+Miro+Monet+Renoir=Artists")
78436+7862553+7219+79436+134921=8162565
> smm("apple+grape+plum=banana")
67723+89673+7250=164646
> smm("seven+seven+six=twenty")
68782+68782+650=138214
> smm("earth+air+fire+water=nature")
67432+704+8046+97364=173546
> smm("five+five+nine+eleven=thirty")
4027+4027+5057+797275=810386
> smm("Saturn+Uranus+Neptune+Pluto=Planets")
127503+502351+3947539+46578=4623971
> smm("five+seven+eleven+twelve+fifteen+twenty=seventy")
3209+59094+969094+819609+3238994+819487=5909487

最新の20件