## prereq:
## 1) How to define a function

## Security level
## Warning: Does not stop execution
## Error: Stop execution

## example of warning

log(-1)

## example of error

log("ABC")

## example of non-warning and non-error, but potentially problematic

log(NA)
mean(c(1,2,43,5,6,7,NA))
10/0

### Real life bugs are not only about error / warning.

## how to throw an error / warning
stop(2 > 3)
stopifnot(3==4)

warning("You should go to the Hong Kong Python usergroup!")

suppressWarnings(log(-1))

## let's define a simple function

simpleFx <- function(x) {
    if (x > 0) {
        return("Positive")
    } else {
        return("Negative")
    }
}

# let's try

simpleFx(10)
simpleFx(-10)
simpleFx(50/10)
simpleFx(pi)
simpleFx(10/0)

# let's try some corner cases

simpleFx(0)

simpleFx(NA)
simpleFx(log(-1))
simpleFx(NaN)

NaN > 0

# let try something even crazy
# press C-c C-c to stop

# finding square root using bisection method

mySqrt <- function(x, epsilon = 0.001) {
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs((guess * guess) - x) > epsilon) {
        if ((guess * guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
    }
    return(guess)
}

### infinite loop! Why?
### oldskool way: adding print statement


mySqrtdebug <- function(x, epsilon = 0.001) {
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs((guess * guess) - x) > epsilon) {
        cat(paste0("guess: ", str(guess)))
        cat(paste0("Highest: ", str(highest)))
        cat(paste0("Lowest: ", str(lowest)))
        cat((paste0("Diff: ", str(abs((guess*guess) - x)))))
        if ((guess * guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
    }
    return(guess)
}


### guess is not updated!

mySqrtdebug <- function(x, epsilon = 0.001) {
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs((guess * guess) - x) > epsilon) {
        print(guess)
#        cat(paste0("Highest: ", str(highest)))
#        cat(paste0("Lowest: ", str(lowest)))
#        cat((paste0("Diff: ", str(abs((guess*guess) - x)))))
        if ((guess * guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
        guess <- mean(c(lowest, highest))
    }
    return(guess)
}


# how to detect the problem without inserting print / cat statement

debug(mySqrt)
mySqrt(5) # c (or just enter), Q, help | print out the value of guess after each iteration | ls()
undebug(mySqrt)
mySqrt(5)


mySqrt <- function(x, epsilon = 0.001) {
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs((guess * guess) - x) > epsilon) {
        if ((guess * guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
        guess <- mean(c(lowest, highest))
    }
    return(guess)
}

debug(mySqrt)
mySqrt(5)
undebug(mySqrt)

mySqrt(6)

mySqrt(100)
mySqrt(25.8) * mySqrt(25.8)

### corner cases

mySqrt(0.000001) # Accuracy problem
mySqrt(-1) # also, it should return a complex number
mySqrt(99999999)
mySqrt(Inf) #die
mySqrt(NaN) #die too
mySqrt(NA) #die three

debug(mySqrt)
mySqrt(-1)
mySqrt(0.000001) # Accuracy problem

# defensive programming

mySqrt <- function(x, epsilon = 0.001) {
    if (!(is.numeric(x) & !is.infinite(x) & x > 0)) stop("Invalid input")
    if (x < epsilon) stop("x is smaller than tolerance.")
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs((guess * guess) - x) > epsilon) {
        if ((guess * guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
        guess <- mean(c(lowest, highest))
    }
    return(guess)
}

debug(mySqrt)
mySqrt(Inf)
mySqrt(NaN)
mySqrt(NA)
mySqrt(-1)
mySqrt(0.000001)

## Example from the SICP

## Newton-Rhapson method

sqrtNewt <- function(x, epsilon=0.0001) {
    if (!(is.numeric(x) & !is.infinite(x) & x > 0)) stop("Invalid input")
    if (x < epsilon) stop("x is smaller than tolerance.")
    guess <- max(c(x,1)) / 2
    delta <- abs(guess^2 - x)
    while (delta > epsilon) {
        guess <- mean(c(guess, (x/guess)))
        delta <- abs(guess^2 - x)
    }
    return(guess)
}

sqrtNewt(Inf)
sqrtNewt(NaN)
sqrtNewt(NA)
sqrtNewt(-1)
sqrtNewt(0.000001)
sqrtNewt(100)
sqrtNewt(500)


newton <- function(x, f, epsilon=0.0001) {
    if (!(is.numeric(x) & !is.infinite(x) & x > 0)) stop("Invalid input")
    if (x < epsilon) stop("x is smaller than tolerance.")
    guess <- max(c(x,1)) / 2
    delta <- abs(f(guess) - x)
    while (delta > epsilon) {
        guess <- mean(c(guess, (x/guess)))
        delta <- abs(f(guess) - x)
    }
    return(guess)
}

sq <- function(x) {
    return(x^2)
}

newton(100, sq)

# or even, aka lambda function

newton(100, function(x) x^2)

newton(100, function(x) x^3)

debug(newton)
newton(100, function(x) x^3)

# the theory is wrong, because the improvement function depends on f(x)
# i.e. Universial guess <- mean(c(guess, (x/guess))) as guess improvement is wrong


newton <- function(x, f, imp, epsilon=0.0001) {
    if (!(is.numeric(x) & !is.infinite(x) & x > 0)) stop("Invalid input")
    if (x < epsilon) stop("x is smaller than tolerance.")
    guess <- max(c(x,1)) / 2
    delta <- abs(f(guess) - x)
    while (delta > epsilon) {
        guess <- imp(x, guess)
        delta <- abs(f(guess) - x)
    }
    return(guess)
}

newton(100, function(x) x^2, function(x, y) { mean(c(y, (x/y))) })

newton(90, function(x) x^3, function(x, y) { ((x/(y^2)) + (2 * y)) / 3 })


bisec <- function(x, f, epsilon = 0.001) {
    if (!(is.numeric(x) & !is.infinite(x) & x > 0)) stop("Invalid input")
    if (x < epsilon) stop("x is smaller than tolerance.")
    lowest <- 0
    highest <- x
    guess <- mean(c(lowest, highest))
    while(abs(f(guess) - x) > epsilon) {
        if (f(guess) > x) {
            highest <- guess
        } else {
            lowest <- guess
        }
        guess <- mean(c(lowest, highest))
    }
    return(guess)
}


bisec(100, function(x) x^2)
bisec(1000, function(x) x^3)
bisec(10000, function(x) x^4)