Write a function, called num_steps
that takes one input called N
, representing the number of discs.
The function must return the number of steps required to move N
discs from peg 𝐴 to peg 𝐵 following the game’s rules.
How do you know this is a valid answer?
Can you give a mathematical proof of your answer?
You can only give this answer at Question 4
If N=0 number of steps=0
If N=2 number of steps=3
If N=3 number of steps=7
If N=4 number of steps=15
If N=5 number of steps=31
If N=6 number of steps=63
Formula of number of steps that depend on N is “(2^N)-1”
print
is not return
Do not print
. Don’t tell secrets. Don’t break the contract.
If you don’t understand, it’s probably my fault
If you don’t ask, then it’s certainly your fault
Can you believe this?
Formula for Fibonacci numbers. If \(F_n=F_{n1}+F_{n-2}\) and \(F_1=F_2=1,\) then \[F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right)\]
How do you know if this is true? Can you test it for \(F_{10^{100}}\)
Let’s say we have \[1, 4, 8, 48, 88, 488, …\]
What comes next? Why?
Now you can apply the num_steps
function to the values of N
from 1 to 30, and store the results in a vector named S
.
You can plot S
depending on N
and you will see this diagram.
This plot will help you to decide which transformation to use for fitting a linear model. Please build a good linear model and print the resulting coefficients.
log(S)
and print the coefficients
log
and exp
Write a function, called hanoi()
, taking four inputs:
N
: the number of discs to moveorigin
: the name of the peg where the discs are nowdestination
: the name of the peg where the discs must be at the endtemporary
: the name of the peg that can be used as a temporary place for some discsThe function must return a vector of text. Each element of the vector has the name of two pegs, indicating the movement of a single disc. The easiest way to make such text is using the paste()
function.
hanoi <- function(N, origin, destination, temporary) {
if (N == 1) {
return(paste(origin, destination))
} else {
return(c(paste(hanoi(N-1, origin, temporary, destination)),
paste(hanoi(1, origin, destination, temporary)),
paste(hanoi(N-1, temporary, destination, origin))))
}
}
We don’t need paste(x)
. Just x
(but paste(x, y)
is useful)
hanoi <- function( N, origin, destination, temporary) {
if( N == 1) {
print(paste(origin, destination))
}
else{
hanoi(N-1, origin, temporary, destination )
print(paste(origin, destination ) )
hanoi(N-1, temporary, destination, origin)}
}
Do not print
. Don’t tell secrets. Don’t break the contract.
hanoi <- function(N, origin, destination, temporary) {
toh_in <- function(n, s, d) {
if (n == 1) {
x[i] <<- paste(temp[1:3 == s], temp[1:3 == d])
i <<- i + 1
}
else {
toh_in(n - 1, s, 6 - s - d)
x[i] <<- paste(temp[1:3 == s], temp[1:3 == d])
i <<- i + 1
toh_in(n - 1, 6 - s - d, d)
}
return(x[(i - 2 ^ n + 1):(i - 1)])
}
temp <<- c("A", "B", "C")
s = which(temp == origin)
d = which(temp == destination)
x <<- i <<- 1
toh_in(n = N, s, d)
}
hanoi <- function(N, origin, destination, temporary) {
origin <- rep(NA, N)
destination <- rep(NA, N)
temporary <- rep(NA, N)
origin[1] <- temporary
destination[1] <- 0
for(i in 2:N) {
origin[i] <- destination[i-1]*temporary
destination[i] <- destination [i-1] + origin[i]
}
return(data.frame(origin, destination))
}
paste(origin, destination)
hanoi <- function(N, origin, destination, temporary) {
origin <- d_origin <- rep(NA, N)
destination <- d_destination <- rep(NA, N)
temporary <- d_temporary <- rep(NA, N)
origin[1] <- N
d_origin <- 0
destination[1] <- d_destination <- 0
temporary[1] <- d_temporary <- 0
for(i in 2:N) {
d_origin[i] <- -temporary[i-1] * destination[N]
d_destination[i] <- + temporary[i-1] * origin[N]
d_temporary[i] <- + origin[i-1] * - destination[i-1]
origin <- origin[i-1] + d_origin[i]
destination <- destination[i-1] + d_destination[i]
temporary <- temporary[i-1] + d_temporary[i]
}
paste(origin, destination, temporary)
}
What is the formula for S
depending on N
? To get an exact result, you may use first a linear model to find the relation between S+1
and N
. How long would it take to move 64 discs if you move 1 every second?
# write your answer here in plain English
When we add 1 to S, our coefficient numbers become 1 and 2.
In every stage, we multiply the previous stage result, and add 1 extra step. So the formula should be something like this: 2 * hanoi(N-1) + 1.
For a 64 discs game, we need to move:
1.844674e+19 times. If we move 1 disc in every second, it will take 1.844674e+19 second. A year, basically has 3e+06 seconds.
It would take 6e+12 year, or 60 billion centuries.
Finding Formula: 1. Plotting S
depending on N
. - We can use plot() function. 2. Building a good linear model. - We should use log() function on S
and use lm() function. We should store the results. 3. Finding coefficients of the good linear model. - We should use exp() function for the stored result. We will get 2 values. 4. We should write y= First value * second value^N. - If we write to all of codes correct, we will get this formula: y=0.8536174*2.0151370^N
Moving 64 Discs:
1-We will find out how many moves we have to move 64 disks. -For this, we can use the function we wrote in the first question. We will get the result 1.844674e + 19. This number gives us the total number of discs we will move. 2-We multiply this number by 1 second.
Finding Formula:
S
depending on N
.
S
and use lm() function. We should store the results.Moving 64 Discs:
I compared the two linear models with the two following codes using the plot:
lines((exp(predict(lm(log(S)~N)))~N))
lines((exp(predict(lm(log(S+1)~N)))~N))
When I compared both outputs, I concluded that the second linear model yields much more accurate results.
Thus, I calculated the coefficients of the second linear model, which are 1 and 2, respectively. Therefore, the formula for S+1
depending on N
is:
S+1=1*2^N
S=2^N-1
It would take 2^64-1 (1.844674e+19) seconds to move 64 discs.
For Every N, it can be calculated as “2 ^ N - 1”. So S is calculated. To move 64 disks, we need to move the disks 1,845 e+19 times. This process takes more than 500 million years.
Formula could be written Tn = 2^n-1.
n = number of discs.
So;
T64 = 2^64-1
= 1.85x10^19
Formulation of Hanoi moves S=(N^2)-1. We can explain if we want to find the number of moving(N) for depends on S, we calculate square root of S+1 for find to N value or we can improve a code for count move discs timing.
find_time<-function(S) {
second<-S%%60
S<-S%/%60
minute<-S%%60
S<-S%/%60
hour<-S%%24
S<-S%/%24
day<-S%%30
S<-S%/%30
month<-S%%12
S<-S%/%12
year<-S
print(paste("if every move take a second it will take ",year,"years ",month,"month", day, "days", hour, "hours",minute, "minutes", second,"seconds"))
}
find_time(num_steps(64))
remove(find_time)
According to the legend, at the beginning of time, the Hindu temple priests were given a tower of 64 disks of gold and they had the task of transferring the disks from one pole to the third pole on the other side of the temple
The priests worked day and night to complete the task, and when they did, legend has it that the temple crumbled into dust, and the world vanished before the priests could complete the task
The legend was spread by Edouard Lucas in hopes of helping to popularize the Tower of Hanoi game. He invented the “Tower of Hanoi” puzzle in 1883
Read E.coli genome (accession NC_000913) and store it in genome
.
To verify, print the genome length
library(seqinr)
read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)
genome <- read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)
Function read.fasta
gives a list
We want the first list element
library(seqinr)
genes <- read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)
genome <- genes[[1]]
Here we are reading genes ("FASTA features nucleotides, .ffn
), not chromosomes (Fasta nucleotides, .fna
)
library(seqinr)
chromosomes<- read.fasta("C:/Users/Hp2020/Downloads/NC_000913.fna", seqtype= "DNA", set.attributes = FALSE)
genome <- chromosomes[[1]]
This does not work in other people’s computers
install.packages("seqinr")
library(seqinr)
setwd("C:/Users/Casper/Desktop")
chromosome <- read.fasta("sequence.fasta" , seqtype = "DNA")
genome <- chromosome[[1]]
Never use install.packages
or setwd
in a Rmarkdown file
The score of a DNA fragment is calculated as the sum of the value of the matrix on for each letter on each position. That is why this kind of matrices are called “Position Specific Scoring Matrix”, or PSSM.
The first letter takes the value from the first row, the second letter looks into the second row, and so on. Then all the values are summed together, that is the score. For example, the sequence ACTGACTGA is \[-186 + -127 + -186 + -286 + -86 + 184 + -127 + -186 + 167\]
Please write an R function, called calc_score
to evaluate the score of a DNA fragment. The function takes two inputs: the DNA fragment, and a Position Specific Scoring Matrix, and returns the score of that fragment. That is, the sum of all values in the corresponding row and column.
calc_score <- function(dna, M) {
V <- toupper(dna)
scores <- rep(NA, length(V))
for(i in 1:length(V)) {
if(V[i] == "A") {
scores[i] <- M$A[i]
} else if(V[i] == "T") {
scores[i] <- M$T[i]
} else if(V[i] == "G") {
scores[i] <- M$G[i]
} else if(V[i] == "C") {
scores[i] <- M$C[i]
}
}
ans <- sum(scores)
return(ans)
}
Using the calc_score
function from the previous question, we will now calculate the score for any position in the genome. Please write a function called score_of_position
that takes three inputs: the position in the genome (a number), the PSSM matrix, and the genome.
score_of_position <- function(i, M, genome) {
Sum <-0
position <- i
genome <- toupper(genome)
for (j in 1:9) {
if (genome[position +j-1] == "A"){
Sum <- Sum+M[j,1]}
if (genome[position +j-1] == "C"){
Sum <- Sum+M[j,2]}
if (genome[position +j-1] == "G"){
Sum <- Sum+M[j,3]}
if (genome [position +j-1] == "T"){
Sum <- Sum+M[j,4]}}
return(Sum)}
score_of_position <- function(i, M, genome) {
x <- toupper(genome[seq(from=i,length.out=9)])
score_vector <- rep(NA, length(x))
for(i in 1:length(x)) {
if (x[i] == "T") {
score_vector[i] <- M$T[i]
} else if(x[i] == "G") {
score_vector[i] <- M$G[i]
} else if(x[i] == "C") {
score_vector[i] <- M$C[i]
} else if(x[i] == "A") {
score_vector[i] <- M$A[i]
}
}
ans <- sum(score_vector)
return(ans)
}
score_of_position <- function(i, M, genome) {
position <- seq(from=i, length.out=9)
window <- genome[position]
answer <- calc_score(window, M)
return(answer)
}
Why 9?
Why 8?
Now it works with any matrix M
In class 9 we used the GC skew to locate a genomic region where the replication origin may be located. The binding site for DnaA should be in the same range of positions.
Now you can apply the score_of_position
function to evaluate the score for each element of position
, and store the resulting values in a new vector called score
. This answer is just code, not a function.
Finally, to find the binding site for DnaA, we need to find the place that has the largest score. We do not care about the score, only about its position in the genome.
To verify that you read this questions at the right time, please send me an email to the answer address with your student number in the subject and the text “I read it”.
Please complete the code of the food_mice_cats
function. The basic parts of the code have already been filled for you. Do not change the code outside the answer region.
birth <- birth_rate * fat_mice[i-1]
growth <- growth_rate * food[i-1]
reprod <- reprod_rate * fat_cats[i-1]
catch <- catch_rate * fat_mice[i-1] * cats[i-1]
death <- death_rate * fat_cats[i-1]
eating <- eating_rate * food[i-1] * mice[i-1]
hunger <- hunger_rate * mice[i-1]
dying <- dying_rate * cats[i-1]
d_food[i] <- growth - eating
d_mice[i] <- birth - eating - hunger
d_cats[i] <- reprod - catch - dying
d_fat_cats[i] <- catch - death
d_fat_mice[i] <- eating - catch
In the previous question we saw that mice eat the food until the cats catch most of the fat_mice. One strategy that mice can use to survive is to eat faster so they can grow faster and outcompete the cats. To explore this question you will write a function, called final_state
, that takes one input called rate
.
This function should simulate the food_mice_cats
system for 180 days using the same initial conditions and rates as in the previous question, except that the eating_rate
input of food_mice_cats
should be made equal to rate
.
The output of food_mice_cats
should be a vector with four elements. The element names should be growth
, food
, total_cats
, and total_mice
. The value of growth
should be equal to rate
. The value of food
must be the quantity of food on the last day of simulation (day 180). The value of total_cats
must be the sum of the number of cats
and fat_cats
on the last day. In the same way, total_mice
must be the sum of the number of mice
and fat_mice
on the last day.
final_state <- function(rate) {
result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
growth <- rate
food <- result$food[180]
total_cats <- result$cats[180] + result$fat_cats[180]
total_mice <- result$mice[180] + result$fat_mice[180]
return(c(growth, food, total_cats, total_mice))
}
final_state <- function(rate) {
result <-
food_mice_cats(
N = 180, birth_rate = 0.1, growth_rate = 1e-3,
reprod_rate = 0.04, catch_rate = 0.05, death_rate = 0.01,
eating_rate = rate, hunger_rate = 1e-3, dying_rate = 0.05,
food_ini = 100, mice_ini = 10, fat_mice_ini = 0,
cats_ini = 10, fat_cats_ini = 0)
result <- tail(result, 1)[c("food", "cats", "fat_cats",
"mice", "fat_mice")]
result <- c(rate, result[1, 1], sum(result[2:3]), sum(result[4:5]))
names(result) <- c("growth", "food" , "total_cats" , "total_mice")
result
}
final_state <- function(rate) {
food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3, reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0,
cats_ini=10, fat_cats_ini=0) -> a
answer_final_state<- c(growth = rate, food = a$food[180],
total_cats = sum(a$cats[180] + a$fat_cats[180]),
total_mice = sum(a$mice[180] + a$fat_mice[180]))
return(answer_final_state)}
final_state <- function(rate) {
result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
return(c(growth = rate,
food = result$food[180],
total_cats = result$cats[180]+result$fat_cats[180],
total_mice = result$mice[180]+result$fat_mice[180]))
}
final_state <- function(rate) {
keeper <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
growth<-rate
food<-keeper["food"][[1]][180]
total_cats<-keeper["fat_cats"][[1]][180]+keeper["cats"][[1]][180]
total_mice<-keeper["fat_mice"][[1]][180]+keeper["mice"][[1]][180]
return(tibble(growth,food,total_cats,total_mice))
}
final_state <- function(rate) {
for(i in 1:length(rate)){
data <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate[i],
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)}
out <- c(rate[i], data[180,2,i], data[180,4,i]+data[180,5,i],
data[180,3,i]+data[180,6,i])
return(out)
}
final_state <- function(rate) {
n<-180
newResult <- food_mice_cats(N=n, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
namer<-list(rate,newResult$food[n],
(newResult$cats[n] + newResult$fat_cats[n]),
(newResult$mice[n] + newResult$fat_mice[n]))
names(namer[[1]])<-"growth"
names(namer[[2]])<-"food"
names(namer[[3]])<-"total_cats"
names(namer[[4]])<-"total_mice"
return(c(namer[[1]],namer[[2]],namer[[3]],namer[[4]]))
}
final_state <- function(rate) {
result2 <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
answer <- c( growth = rate, food = result2$food[180], total_cats = sum(result2$cats[180], result2$fat_cats[180]), total_mice = sum(result2$mice[180], result2$fat_mice[180])
)
answer
}
final_state <- function(rate) {
var1 <-tail(
food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3, reprod_rate=0.04,
catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0
), 1)
food <- var1$food
total_cats <- var1$cats + var1$fat_cats
total_mice <- var1$mice + var1$fat_mice
return (as.double(data.frame(rate, food, total_cats, total_mice)))
}
final_state <- function(rate) {
temporary_result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
reprod_rate=0.04, catch_rate=0.05, death_rate=0.01, eating_rate=rate,
hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
mice_ini=10, fat_mice_ini=0, cats_ini=10,fat_cats_ini=0)
growth<-rate
food<-temporary_result$food[180]
total_cats<-temporary_result$cats[180] + temporary_result$fat_cats[180]
total_mice<-temporary_result$mice[180] + temporary_result$fat_mice[180]
names(growth)<-"growth"
names(food)<-"food"
names(total_mice)<-"total_mice"
names(total_cats)<-"total_cats"
return(c(growth,food,total_cats,total_mice))
}