I have always loved symmetry in photography. The picture above was created in Photoshop simply by coping and rotating certain parts. See SchlockySymmetry for other symmetrical photos. Recently, I stumbled upon an article in R-bloggers that used abind and biOps for some simple image manipulation. I thought it would be fun to perform some of the tricks mentioned in the article to achieve the picture above.

R Code:

#Image processing

#Read Original Image
x <- readJpeg("Boudha.jpg")

plot(x)
dev.copy(png,filename="original.png",width=2048,height=1206);
dev.off ();

#copy pixels from left to right and mirror
y <- imagedata(abind(x[,1:(ncol(x)/2),], x[,(ncol(x)/2):1,] , along=2))

plot(y)
dev.copy(png,filename="halfmirror.png",width=2048,height=1206);
dev.off ();

#copy all pixels and mirror rows
z <- imagedata(abind(y[1:(nrow(x)),,], y[(nrow(x)):1,,] , along=1))

plot(z)
dev.copy(png,filename="final.png",width=2048,height=1206);
dev.off ();

Result

original.png
halfmirror.png
final.png
"Little boxes on the hillside,
Little boxes made of ticky tacky,
Little boxes on the hillside,
Little boxes all the same.
There's a green one and a pink one 
And a blue one and a yellow one,
And they're all made out of ticky tacky
And they all look just the same."














The Zorganian Republic has some very strange customs. Couples only wish to have female children as only females can inherit the family's wealth, so if they have a male child they keep having more children until they have a girl. If they have a girl, they stop having children. What is the ratio of girls to boys in Zorgania?

The ratio of girls to boys in Zorgania is 1:1. This might be a little counter-intuitive at first. Here are some ways of tackling this problem.

1. Monte Carlo Simulation:

Although, Monte Carlo simulation does not necessarily show why the result is 1:1, it is appropriate because of the very counter-intuitive nature of the problem. At the very least, it helps us see that the result is indeed 1:1. Therefore, this is a good start.

The following R code estimates the probability of a child being a boy in Zorgania. 

couples <- 100000
boycount <- 0

for (i in 1:couples){
    # 0: boy
    while (sample(c(0,1),1) == 0) {
    boycount=boycount+1
  }
}
probability <- boycount/(couples+boycount)

Result:


2. Understanding the question better:

Here is another question: What is the probability of getting a tail in a fair coin toss, if all seven previous tosses resulted in heads? Since coin flips are independent events, the probability is still going to be 0.5. Similarly in this case, the child births are independent. It does not matter if the couples stop giving birth after they have a baby-girl. The expected value is unchanged. 

For example, consider five couples: C1, C2, C3, C4 and C5.
If B-> Boy and G-> Girl. Using R's sample().
For C1:  {B, G}
For C2: {G}
For C3: {B, B, G}
For C4: {B, G}
For C5: {G}

Now, ignore the couples and only consider the children. The children are {B, G, G, B, B, G, B, G, G}. The only thing happening here is simply the generation of  either a B or a G with equal probability for each generation. At this point, it is quite obvious that the part that has to do with "couple...." in the question is to mislead and confuse similar to the "previous seven tosses.." example that I mentioned in the beginning of 2. 

3.  Counting:

If we start with 512 couples (hence 512 first borns), half of them are going to have a girl as their first. Those couples will stop having children. Among, the other half couples who had a son as their first child, half of them are going to have a girl as their second child and so on. At every step there is an equal numbers of boys and girls. Therefore, the expected ratio is 1:1.



This is the naive/brute-force implementation of the Mandelbrot Set plotting. I just followed the algorithm.  

# Plotting the Mandelbrot Set 

# length of sequence for real and imaginary parts of complex numbers
length <- 1000

# sequences for real and imaginary parts
real = seq(-1.8,0.6, len=length)
imaginary = seq(-1.2,1.2, len=length)


result <- matrix(nrow = length, ncol = length)

for (i in 1:length)
{
  for (j in 1:length)
  {
    result[i,j]=inmandelbrotset(complex(real = real[i], imaginary = imaginary[j]))       
  }
}

image(result, axes=FALSE)


# function that checks if a point E mandelbrot set
inmandelbrotset <- function(c)
{
  dwell.limit <- 2048
  z <- 0  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}

Adding colors:

We now have a Boolean matrix that records if a point is in the Mandelbrot Set. Since the matrix can only have two values : true or false, thus far, we have only been able to plot read and white images. The next step is to add colors such that we get more information on when a particular point escapes the radius of 2. Again, this is the naive/brute force way of doing it.

# Mandelbrot Plotting with colors 
length <- 1000
real = seq(-2.0,2.0, len=length)
imaginary = seq(-2.0,2.0, len=length)
result <- matrix(nrow = length, ncol = length)
dwell.limit <- 512

for (i in 1:length)
{
  for (j in 1:length)
  {  
    z <- 0
    c <-complex(real = real[i], imaginary = imaginary[j])
    for (k in 1:dwell.limit)
    { 
      
      z <- z ** 2 + c
      if (Mod(z) > 2)
      {      
        result[i,j]=k
        break
      }
    }     
    
  }
}
set.seed(2)
image(result,breaks=0:dwell.limit
      ,col=c(1,sample(terrain.colors
                      (dwell.limit-1,alpha = .8))),asp=1,ax=F)







and, just for the heck of it..

ASCII Mandelbrot Set using R (naive)


s <- seq(-1.7,1.2, by =.1)
a <- ""
for (i in 1:length(s))
{  
  for (j in 1:length(s))
  {
   a<-cat(a,inmandelbrotset(complex(r = s[j], i = s[i])))     
  }  
  a <- cat(a,"\n")
}





Achieved by returning a " " or "#" instead of FALSE or TRUE from function "inmandelbrotset".

A better algorithm

Utilizing R's easy to use lists in implementation:



# more efficient algorithm to plot the Mandelbrot set 

sequence <- seq(-2,2,len=1000)
dwell.limit <- 200

# matrix of points to be iterated 
complex.matrix <- t((sapply(sequence,function(x)x+1i*sequence)))
in.mandelbrot.index <- 1:length(complex.matrix)
iter=z=array(0,dim(complex.matrix)) 
           
for(i in 1:dwell.limit){ 
  # complex quadratic polynomial function for all points
  z[in.mandelbrot.index]=complex.matrix[in.mandelbrot.index]+z[in.mandelbrot.index]^2 
  # boolean matrix
  result=Mod(z[in.mandelbrot.index])<=2  
  # if result is false, store the iteration 
  iter[in.mandelbrot.index[!result]]=i
  # save all the index where points are still in the mandelbrot
  in.mandelbrot.index=in.mandelbrot.index[result]
}
set.seed(19)
image(iter,main=paste("Iterations: ", i, sep=" "), breaks=0:dwell.limit
      ,col=c(1,sample(rainbow
                      (dwell.limit-1,alpha = .8))),ax=F, asp=1)

Plotting the Julia set

A little modification to the code above (red and white Mandelbrot) produces Julia Sets. The idea here is to set a constant C and send Z to the function instead of C.

c <- complex(real=-0.1,imaginary=0.651) 
label <- toString(c)
injulia <- function(z)
{
  dwell.limit <- 128
  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}



Adding colors:

This is achieved by following the same process as above.






Sierpinski Gasket using Chaos game

#### Chaos game for generation of Sierpinski Gasket
# 1. Take 3 points in a plane to form a triangle, you need not draw it.
# 2. Randomly select any point inside the triangle and consider that your current position.
# 3. Randomly select any one of the 3 vertex points.
# 4. Move half the distance from your current position to the selected vertex.
# 5. Plot the current position.
# 6. Repeat from step 3

plot.new()
iterations <- 2000
vertices <- matrix(c(0,0,0.5,1,1,0),3,2, byrow=T)
current.point <- c(0.5,0.5)
random.vertex <- sample(1:3,iterations,replace=T)
plot.result = matrix(nrow=iterations,ncol=2)
for (i in 1:iterations){
 current.point <- (current.point+vertices[random.vertex[i],])/2
 plot.result[i,] <- current.point
}
points(plot.result,pch = 46)



Adding colors:
points(plot.result,pch = 46,col=c(13,3,41)[random.vertex])


I used R for two of my classes in college: Probability and Statistics and Stochastic Probabilities. During that time, I was in reverence with the way it handled lists. However, I did not use it much after getting done with those classes. I have decided to relearn it.

In this post, I will start off with some simple programs. After ascertaining some comfort with the syntax and semantics of R, I will try some complex ones. Also, since I am in the process of learning myself, if you decide to use some of these codes, just keep in mind that there might be several better ways of doing the same thing.

1. Dice Rolls: Probability of getting a side f number of times in n dice rolls

# DiceRolls.R
# Probability of getting a side f number of times in n dice rolls

#inputs
side <- 4
f <- 2
n <- 5

# Binomial Distribution
# using R's function for binomial distribution 
# where P(X=k)
 probability <- dbinom(frequency, n, 1/6)
 probability

# Monte Carlo Method
rep <- 10000 # Number of repetitions
list <- numeric(rep) # Initialize vector
for (i in 1:rep) {
trial <- sample(1:6, n, replace=TRUE)
#1 if number of times a side appears is f
success <- if ((sum(trial==side))==f) 1 else 0
list[i] <- success 
}
mean(list)

>  probability
[1] 0.160751
> mean(list)
[1] 0.1647

2. One dimentional random walks

# Five random walks in one dimension starting at 0

# creating a color string vector
color <- c("red","green","blue","black","darkgoldenrod1")
# cleaning the plot; this is necessary because we are using par(new=TRUE)
plot.new()
for (i in 1:5) {
x0<-0
I<-1
T <-1000
x=c(x0,rnorm(T-1)) # T-1 randomly generated normal distributions
y=cumsum(x)
par(new=TRUE) # so that the previous graphs are not deleted
plot(y,type="l",col=color[i], main="Five random walks in one dimension starting at 0", xlab="", ylim=c(-60, 60))
}




3. Two dimensional random walk


I stumbled upon a presentation slideshow that discusses how R can be used more efficiently. They perform analysis on several implementations that simulate two dimensional random walk. I started out with what they call "naive implementation". According to them, this is how a person who codes in C, C++ or Java would think of doing it. I totally agree. This exercise is great because it helps understand the positives of R in comparison to other languages that I am comfortable with. Link to the aforementioned presentation slideshow. 

Naive implementation:
# naive implementation (loops are slow in R)
naive <- function(n){

# numeric creates a double-precision vector of the specified length with each element equal to 0
xpos = ypos = numeric(n)

# will be used to randomly generate direction
xdir = c(TRUE,FALSE)
nextdir = c(1,-1)

# xpos and ypos are already 0 (see comment on numeric)
for(i in 2:n){

 if(sample(xdir,1)){
  xpos[i]=xpos[i-1]+sample(nextdir,1)
  ypos[i]=ypos[i-1]
 }
 else{
  xpos[i]=xpos[i-1]
  ypos[i]=ypos[i-1]+sample(nextdir,1)
 }
}
plot(xpos,ypos,type="l",col="blue",main="Two Dimensional Random Walk", xlab="x", ylab="y")
# starting point
points(0,0,pch = 19,col="red")
}

> system.time(naive(50000))
   user  system elapsed 
   0.84    0.17    1.02 


> system.time(naive(100000))
   user  system elapsed 
   1.70    0.32    2.03 
Vector implementation:

This implementation gets rid of the for loop. It is similar to the code in one dimensional random walk code above.

# Vectorization Process
factorization = function(n){ 
 steps = sample(c(-1,1), n-1, replace = T)
 xdir=sample(c(TRUE,FALSE), n-1,replace=T)
 xpos=c(0,cumsum(ifelse(xdir,steps,0)))
 ypos=c(0,cumsum(ifelse(xdir,0,steps)))
 plot(xpos,ypos,col="green",type="l",main="Two Dimensional Random Walk", xlab="x", ylab="y")
 points(0,0,pch = 19,col="red")
}

> system.time(vectorization(50000))
   user  system elapsed 
   0.05    0.19    0.25 


> system.time(vectorization(100000))
   user  system elapsed 
   0.09    0.31    0.40 

4. Birthday Probability

Click here to go to my post on Birthday Probability using Java. Using R, I was really surprised how easy it was.

How a java programmer would have done it in R

# Birthday Probability
birthday <- function(numofpeople, numofrooms){
  count <- 0
  for (i in 1:numofrooms){
    bday <- sample(1:365,numofpeople, replace=T)
    if (any(duplicated(bday)))
      count<-count+1;
  }  
probability = count/numofrooms
probability
}

How a R programmer would have done it (I think)

# Birthday Probability
birthday <- function(numofpeople, numofrooms){ 
  list <- numeric(numofrooms)
  
  for (i in 1:numofrooms){
    bday <- sample(1:365,numofpeople, replace=T)
    list[i] <- if (any(duplicated(bday))) 1 else 0   
  }   
mean(list)
}

Result:

> source("birthday.r")
> birthday(23,5000)
[1] 0.5044

5. Mandelbrot Set Area

Click here to go to my previous blog post on the area of the Mandelbrot Set using Java and Javascript. 

monte.Carlo <- 5000
x <- runif(monte.Carlo, -2, 2)
y <- runif(monte.Carlo, -2, 2)
list <- numeric(monte.Carlo)

for (j in 1:monte.Carlo){
  list[j] <- if (inmandelbrotset(complex(real = x[j], imaginary = y[j]))) 1 else 0
}
area<-mean(list)*16

# function that checks if a point E mandelbrot set
inmandelbrotset <- function(c)
{
  dwell.limit <- 2048
  z <- 0  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}

Result:

> system.time(source("mandelbrot.r"))
   user  system elapsed 
   1.91    0.00    1.90 
> area
[1] 1.5168

These are screenshots of some of my very first apps. They did not do much. However, they build my appreciation towards Android development.